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separate  into  isolated  protofilaments . 97 

B.l  A  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely  sep¬ 
arated  chains  of  length  y/N)  of  several  different  patterns  of  h  =  0.001, 
w  =  0.25  particles  versus  y/N .  All  patterns  except  the  ring  pattern 
are  rectangular  lattices  of  dimensions  N  particles  by  N  particles. 

Note  that  the  microvortex  and  macrovortex  energies  are  nearly  iden¬ 
tical  and  are  the  energetically  preferred  patterns  for  y/N  >  3.  Also 
note  that  the  HCP  ferromagnetic  pattern  decreases  significantly  rel¬ 
ative  to  all  other  patterns,  over  the  calculated  range . 108 
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B.2  A  plot  of  the  energies  (normalized  to  the  energy  of  \/rN  infinitely 
separated  chains  of  length  v^A)  of  several  different  patterns  of  h  = 
0.001,w  =  0.31  particles  versus  a /N.  All  patterns  except  the  ring  pat¬ 
tern  are  rectangular  lattices  of  dimensions  N  particles  by  N  particles. 
Note  that  the  microvortex  is  the  energetically  preferred  pattern  for 
a//V  >  3.  Also  note  that  the  HCP  ferromagnetic  pattern  decreases 
significantly  relative  to  all  other  patterns,  over  the  calculated  range. 
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1.  INTRODUCTION 


In  this  thesis,  we  describe  results  from  investigations  of  the  self-assembly  of  anisotrop¬ 
ically  interacting  particles.  In  particular,  we  focus  upon  the  roles  of  dipolar  and 
higher  order  multipolar  interactions  on  the  patterns  of  self-organization.  Using  an 
experimental  model  system  of  vertically  vibrated  magnetic  spheres,  we  investigate 
the  effects  of  octopolar  and  higher  order  interactions  on  the  pattern  of  self-assembly. 
We  show  that  simple  theoretical  point  charge  models  can  be  used  to  provide  insight 
into  the  underlying  causes  of  the  observed  phenomena.  We  also  show  that  such 
models  can  be  used  to  better  understand  the  pattern  formation  in  several  related 
physical  systems,  including  biological  macromolecular  self-assembly  and  cohesive 
granular  materials. 


1.1  Motivation 

The  self-assembly  of  anisotropically  interacting  particles  is  currently  a  topic  of  con¬ 
siderable  interest  for  several  fields  of  science.  In  particular,  the  self-assembly  of 
biological  macromolecules  (which,  in  general,  are  anisotropically  interacting)  is  of 
significant  interest  to  the  pharmaceutical  industry  because  of  possible  therapeutic 
value  [1,  2,  3,  4,  5].  For  example,  the  popular  chemotherapy  drug  Taxol  targets 


the  self-organization  of  tubulin  proteins  in  cancer  cells  and  prevents  cell  division 
and  the  growth  of  tumors  [2,  5].  Understanding  the  conditions  which  are  favorable 
to  self-assembly  and  what  causes  such  conditions  to  be  favorable  is  fundamental  to 
the  control  of  the  polymerization.  As  a  basic  step  toward  understanding  the  self- 
assembly  of  anisotropically  interacting  particles,  we  investigate  the  self-organization 
of  spheres  with  dipolar  and  higher  order  interactions. 

Two-dimensional  dipolar  systems  have  been  investigated  using  many  differ¬ 
ent  experimental  and  theoretical  approaches  [6,  7].  Patterns  have  been  observed 
experimentally  in  nm  to  nun-sized  magnetic  particles  suspended  on  the  surface 
of  a  liquid  [8,  9,  10,  11,  12],  in  a  magnetic  fluid  containing  nonmagnetic  parti¬ 
cles  [10,  13],  in  magnetic  bubbles  in  garnet  films  [14],  in  mm  to  cm-sized  mag¬ 
netically  excited  magnetized  spheres  [15],  at  the  atomic  scale  in  thin  magnetic 
films  [16,  17,  18,  19,  20,  21,  22,  23],  in  dipolar  molecular  fluids  such  as  hydrogen 
fluoride  [24,  25],  and  recently  in  mm-sized  vertically  vibrated  magnetic  spheres  [26]. 
Studies  of  biological  systems  also  suggest  that  dipole-dipole  interactions  between 
biological  molecules  may  be  an  important  part  e.g.  of  pattern  formation  in  a  lipid 
monolayer  [27,  28].  Monte  Carlo  simulations  of  the  thermal  equilibrium  of  large 
2D  systems  of  hard  spheres  with  centrally  embedded  point  dipoles  have  shown  that 
ordering  in  the  absence  of  external  forces  depends  on  the  dipolar  strength,  the  tem¬ 
perature  of  the  system,  and  the  number  density  of  the  spheres  [11,  29,  30,  31]. 

Significant  work  has  also  been  done  on  pattern  formation  in  3D  systems  of 
dipolar  spheres  [32,  33,  34,  35].  The  investigations  show  that  low-density  3D  sys¬ 
tems  of  dipolar  spheres  form  chains  and  rings  at  low  temperatures,  but  the  chains 
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and  rings  do  not  coagulate  to  form  a  liquid  phase.  The  lack  of  a  liquid  phase  was 
suggested  to  be  a  consequence  of  the  saturation  of  the  interaction  potential  during 
chain  formation,  and  the  resulting  lack  of  strong  interactions  between  chains.  In 
addition,  because  of  the  dissociation  and  recombination  of  the  chains  during  equi¬ 
librium,  this  system  of  dipolar  spheres  has  been  compared  to  microscopic  physical 
systems  of  so-called  ‘living  polymers’,  which  exhibit  similar  dynamic  behavior  [32], 
The  focus  of  our  work  is  on  the  self-assembly  of  2D  systems  of  particles  with  dipolar 
and  higher  order  magnetic  interactions.  We  find  that  at  low  temperatures  chains 
can  collect  to  form  a  dense  structure,  in  agreement  with  simulations  of  2D  systems 
of  dipolar  spheres  by  Clarke  et  al.  [25]. 

The  self-assembly  of  anisotropically  interacting  particles  is  of  particular  in¬ 
terest  in  molecular  biology  because  of  the  wide  array  of  biological  systems  which 
exhibit  self-assembly  and  the  variety  of  self-organized  structures  that  are  vital  to 
the  function  of  living  cells  [36,  37].  Studies  of  biological  systems  suggest  that  dipole- 
dipole  interactions  between  particles  play  an  important  role  in  such  phenomena  as 
microtubule  formation  [38]  and  pattern  formation  in  a  lipid  monolayer  [27].  While 
the  interactions  between  these  nanoscale  biological  particles  are  difficult  to  rigor¬ 
ously  model,  simple  model  systems  of  dipolar  spheres  may  provide  insight  into  the 
dynamics  of  biological  self-assembly  in  a  fundamental  manner. 

There  also  is  significant  current  interest  in  the  dynamics  of  granular  materi¬ 
als  that  have  interactions  in  addition  to  the  typical  hard-sphere  interactions.  The 
dynamics  of  wet  granular  mixtures  has  been  a  recent  focus  as  an  example  of  pat¬ 
tern  formation  among  isotropically  attractive  particles  [39,  40].  In  granular  systems, 
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the  addition  of  attractive  forces  between  particles  has  been  observed  to  both  cause 
and  prevent  segregation  under  different  circumstances  [39].  The  dynamics  of  dipo¬ 
lar  interacting  spheres  is  also  of  current  interest  as  a  novel  granular  material  with 
anisotropic  interactions  [41,  42,  26].  Blair  and  Kndrolli  [26]  found  that,  in  a  ver¬ 
tically  vibrated  2D  mixture  of  magnetic  and  non-magnetic  granular  spheres,  the 
magnetic  spheres  can  self  organize  and  cluster,  depending  on  the  volume  fraction 
of  magnetic  particles  and  on  the  relative  strength  of  the  magnetic  dipole-dipole  in¬ 
teraction  to  the  external  vibration  amplitude.  In  our  earlier  work  [41],  we  showed 
that  the  shape  of  the  magnetic  held  is  another  important  factor  in  the  self-assembly, 
ans  the  higher  order  moments  of  the  interactions  can  significantly  affect  the  pattern 
formation. 

Moments  of  order  higher  than  dipole  are  of  particular  importance  to  biolog¬ 
ical  macromolecules,  where  macromolecular  interaction  potentials  can  be  highly 
anisotropic  to  higher  than  dipole  order,  and  small  changes  in  the  interactions  can 
dramatically  alter  the  self-assembly  pattern  (for  example,  see  the  self-assembly  of 
sickle  cell  hemoglobin  in  Chapter  6).  However,  while  the  self-assembly  of  dipolar 
systems  have  been  studied  relatively  extensively  [6,  7,  8,  9,  10,  11,  12,  26,  29,  30, 
31,  34,  43],  the  affects  of  higher  order  moments  on  self-assembly  are,  to  a  large 
degree,  unknown.  Relatively  few  simulations  or  theoretical  studies  using  models 
which  include  higher  order  moments  [24,  25,  27,  54]  have  been  done.  The  purpose  of 
this  thesis  is  a  detailed  study  of  how  higher  order  moments  (in  particular,  octopole 
moments)  affect  the  self-assembly  of  dipolar  particles. 
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1.2  Outline  of  Thesis 


In  Chapter  2,  we  describe  our  experimental  and  theoretical  methods.  We 

give  general  schematics  of  our  experimental  setup,  including  a  general  description  of 
the  dipolar  particles  used.  We  also  describe,  in  general,  the  models  used  to  provide 
theoretical  insight  into  the  observed  self-assembly. 

In  Chapter  3,  we  investigate  the  clustering  transition  as  a  function  of  tem¬ 
perature.  We  show,  using  our  experimental  system  and  Monte  Carlo  simulations 
performed  by  collaborators  at  NIST ,  that  at  temperatures  above  the  clustering  tran¬ 
sition,  the  particles  are  largely  dissociated  from  each  other  and  exhibit  traditional 
gas  phase  characteristics.  At  temperatures  below  the  polymerization  transition,  the 
thermal  (kinetic)  energy  of  the  particles  is  insufficient  to  overcome  the  attractive 
interactions  and  the  particles  form  strongly  bound  clusters.  The  anisotropic  nature 
of  the  dipolar  interactions  induces  head-to-tail  alignment,  resulting  in  the  formation 
of  linear,  ring,  or  branched  structures.  Between  these  two  limits  lies  a  broad  transi¬ 
tion  region  in  which  these  polymer-like  structures  continuously  form  and  disappear 
in  dynamic  equilibrium.  Our  results  show  that  the  transition  temperature,  T$,  in¬ 
creases  and  the  width  of  the  transition  broadens  with  increasing  concentration,  C, 
in  accord  with  the  analytic  theory  of  equilibrium  polymerization.  (Work  to  be  sub¬ 
mitted  to  Physical  Review  E  by  J.Stambaugh,  K.  Van  Workum,  J.  Douglas,  and  W. 
Losert.) 

Chapter  4  describes  the  effects  of  higher  order  multipole  interactions. 

Pattern  formation  is  investigated  for  a  vertically  vibrated  monolayer  of  magnetic 
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spheres.  The  spheres,  of  diameter  D,  encase  cylindrical  magnetic  cores  of  length,  l. 
For  large  D/l,  we  find  that  the  particles  form  a  hexagonal-close-packed  pattern  in 
which  the  particles’  dipole  vectors  assume  a  macroscopic  circulating  vortical  pattern. 
For  smaller  D/l,  the  particles  form  concentric  rings.  The  static  configurational  mag¬ 
netic  energy  (which  depends  on  D/l)  appears  to  be  a  determining  factor  in  pattern 
selection,  even  though  the  experimental  system  is  driven  and  dissipative.  (Work 
published  by  J.Stambaugh,  D.  Lathrop,  E.  Ott,  and  W.  Losert,  Physical  Review  E, 
68  026207  (2003).) 

In  Chapter  5,  we  investigate  binary  systems  of  magnetic  spheres.  We  show, 
using  a  simple  experiment  and  numerical  calculations,  that  in  dilute  self-assembled 
networks,  a  small  change  in  the  multipole  moment  distribution  of  a  small  fraction  of 
the  particles  can  significantly  increase  the  branching  in  the  self-assembled  network. 
We  also  show  that  otherwise  identical  systems  of  magnetic  spheres  can  segregate 
due  only  to  differences  in  their  magnetic  fields.  Vertically  vibrated  monolayers 
of  magnetic  spheres  can  segregate  both  by  field  strength  and  by  field  shape.  In 
binary  systems  of  particles  with  differing  field  strength  and  all  other  properties 
identical,  the  segregation  increases  with  the  proportion  of  weaker  particles,  and 
also  increases  approximately  linearly  with  acceleration  over  the  acceleration  range 
studied.  Binary  systems  of  particles  with  differing  field  shape  and  all  other  properties 
identical  also  show  approximately  linear  increases  in  segregation  as  acceleration  is 
increased  over  the  acceleration  range  studied.  Segregation  occurs  in  conjunction 
with  a  decrease  in  magnetic  energy,  with  the  energy  decrease  being  mostly  due 
to  the  actual  segregation,  rather  than  an  evolution  of  the  spatial  pattern.  (Work 
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published  J.Stambaugh,  Z.  Smith,  E.  Ott,  and  W.  Losert,  Physical  Review  E,  70 
031304  (2004).) 

In  Chapter  6,  we  discuss  biological  self-assembly.  Specific  examples  of  several 
relevant  biological  self-assembling  macromolecnlar  systems  are  given.  We  discuss 
the  important  interactions  which  govern  the  self-assembly,  and,  in  particular,  we 
describe  the  role  of  electrostatic  interactions  in  self-assembly.  We  characterise  the 
electrostatic  properties  of  several  of  these  macromolecules  by  calculating  their  elec¬ 
trostatic  multipole  moments  up  to  octopole.  We  also  describe  a  Monte  Carlo  method 
of  developing  minimal  point  charge  models  of  the  macromolecules  which  match  the 
electrostatic  multipole  moments  of  real  proteins  up  to  the  octopole  terms. 

Chapter  7  describes  a  minimal  interaction  model  of  the  self-assembling 
protein,  tubulin.  We  propose  a  simple  point  charge  model  of  the  protein  tubulin, 
and  we  show  that,  without  altering  the  dipole  moment,  the  addition  of  higher  order 
multipole  moments  to  the  model  can  energetically  stabilise  the  biologically  correct 
self-assembly  pattern.  The  model  consists  of  four  static  point  charges  embedded 
within  an  impenetrable  rectangle  with  the  dimensions  of  a  tubulin  dimer.  We  find 
that  aggregates  of  the  model  tubulin  dimers  reproduce  key  features  of  microtubule 
structure:  cylindrical  shape,  staggering,  and  fraying  at  the  edges.  We  also  present 
a  physical  macroscopic  magnetic  analog  of  our  model,  and  show  that  it  has  many 
of  the  same  structural  characteristics  as  microtubules.  (Work  to  be  submitted  to 
Biophys.  J.  by  J.Stambaugh,  D.  Jones,  and  W.  Losert.) 

Appendix  A  is  a  detailed  calculation  of  the  multipole  moments,  up  to 
octopole,  of  our  model  and  experimental  particles  from  Chapter  5.  This 
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includes  a  determination  of  the  conditions  for  matching  the  moments  of  the  model 
and  experimental  particles. 

In  Appendix  B,  we  theoretically  analyze  the  preferred  self-assembly  pat¬ 
terns  of  some  of  our  magnetic  particles  from  Chapters  4  and  5  as  a  func¬ 
tion  of  the  system  size.  This  includes  a  calculations  of  the  energies  of  several 
different  patterns  of  particles  with  dipolar  and  higher  order  interactions  versus  the 
number  of  particles  in  the  patterns. 

Appendix  C  is  a  detailed  description  of  some  of  the  energy  calculations 
from  Chapter  4.  This  includes  a  calculations  of  the  energies  of  several  different 
patterns  of  particles  with  dipolar  and  higher  order  interactions. 


2.  METHODS 


2. 1  Experimental 

The  following  is  a  description  of  general  characteristics  of  our  experimental  setup. 
Experimental  details  that  are  unique  to  each  chapter  will  be  described  within  the 
respective  Experimental  Details  sections  of  each  chapter.  The  experimental  portion 
of  this  thesis  is  focused  upon  macroscopic  two-dimensional  systems  of  vertically 
vibrated  magnetic  spheres.  Figure  2.1  is  a  schematic  of  the  experimental  setup  used. 
A  horizontal  cylindrical  container  with  a  diameter  17.5  times  a  particle  diameter 
and  a  height  of  1.7  particle  diameters  was  constructed  with  a  black  Delrin  or  gray 
polyvinyl  chloride  bottom  and  a  rigid  transparent  acrylic  top  (for  imaging).  The 
system  was  imaged  from  above  at  a  rate  of  1-2  Hz  using  a  grayscale  or  color  CCD 
camera.  The  acceleration  amplitude,  a,  was  varied  from  4 g  to  11  g  (where  g  is 
the  acceleration  of  gravity)  and  was  measured  using  an  accelerometer.  Magnetic 
spheres  were  placed  in  the  container  and  sinusoidally  vibrated  vertically  at  30  Hz 
using  an  electromagnetic  shaker.  The  maximum  strength  of  the  magnetic  field  from 
the  shaker  at  the  inner  container  surface  was  23.1  G  at  zero  acceleration  amplitude, 
although  at  typical  acceleration  amplitudes  the  RMS  magnetic  field  strength  was  of 
order  1  G.  No  noticeable  effect  due  to  any  external  magnetic  field  from  the  shaker 


was  observed. 


The  magnetic  particles  used  were  cylindrical  permanent  magnets  rigidly  en¬ 
cased  in  spherical  hard  plastic  shells.  The  plastic  shells  and  some  of  the  magnets  (the 
‘long’  magnets)  are  distributed  as  ‘Magnetic  Marbles’  by  the  Safari  toy  company, 
and  were  purchased  in  the  Smithsonian  Air  and  Space  Museum.  Other  magnets 
with  differing  properties  were  purchased  elsewhere  and  were  used  to  replace  the  ex¬ 
isting  Safari  magnets  in  order  to  alter  the  multipolar  interactions  of  the  particles. 
Magnet  materials  were  chosen  for  their  resistance  to  demagnetization  (both  from 
vibration  and  opposing  magnetic  fields)  and  their  physical  strength.  Magnetic  field 
magnitudes  of  our  particles  were  determined  by  measuring  the  maximum  magnetic 
field  on  the  surfaces  of  at  least  twenty  particles  using  a  Gauss  meter  and  calculating 
mean  and  standard  deviation  of  the  measurements.  Figure  2.2  is  a  general  schematic 
of  our  magnetic  particles.  By  altering  the  material  makeup  and  the  axial  length  l  of 
the  magnetic  cores,  we  are  able  to  independently  change  the  dipolar  and  octopolar 
moments  (from  symmetry  the  quadrupole  moments  are  identically  zero).  It  may 
also  be  possible  to  construct  similar  magnetic  particles  with  nonzero  quadrupole 
moments  by,  e.g.,  slicing  a  cylindrical  magnet  vertically  along  its  axis,  rotating  one 
of  the  halves  180°  about  the  axis  perpendicular  to  the  slicing  plane,  then  fixing  the 
halves  back  together.  The  quadrupole  moment  would  then  be  a  function  of  /,  and  the 
particles  would  have  zero  dipole  moment.  Another  method  of  constructing  partlicle 
with  non-zero  quadrupole  would  be  to  permenantly  attach  the  north  faces  of  two 
cylindrical  magnets  together  such  that  the  magnetic  dipoles  are  opposing.  However, 
for  the  remainder  of  this  thesis,  we  will  describe  observations  of  magnetic  particles 
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Fig.  2.1:  A  schematic  of  the  experimental  setup.  Magnetic  spheres  are  placed  in  a  2D 
container  and  vertically  vibrated  at  30  Hz.  The  system  is  imaged  from  above. 


with  zero  quadrupole  moments.  For  the  majority  of  the  experiments  described  in 
this  thesis,  the  particles  were  constructed  to  have  differently  colored  northern  and 
southern  hemispheres  so  that  dipole  orientations  could  be  discerned. 


2.2  Theoretical  Models 

The  following  is  a  description  of  general  characteristics  of  our  theoretical  magnetic 
particle  models.  Model  details  that  are  unique  to  each  chapter  will  be  described 
within  the  respective  Energetic  Considerations  sections  of  each  chapter.  The  theo- 
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Fig.  2.2:  A  general  schematic  of  the  magnetic  particles  used.  The  material  makeup  and 
l  of  the  magnetic  cores  were  parameters  that  were  varied  to  independently  alter 
the  strengths  of  the  dipolar  and  octopolar  interactions. 

retical  models  in  this  thesis  consist  of  two  to  eight  static  point  charges  surrounded 
by  appropriately  shaped  impenetrable  boundaries.  Figure  2.3  is  a  general  schematic 
of  my  model  of  the  magnetic  spheres  described  in  fig.  2.2.  The  model  consists 
of  impenetrable  spheres  of  diameter,  D ,  containing  four  magnetic  ‘charges’  of  equal 
magnitude  ±q,  separated  by  distances  h  and  w  (0  <  h,  w  <  D),  such  that  the  charges 
are  at  the  corners  of  a  rectangle.  This  model  is  a  two-dimensional  theoretical  analog 
of  our  three-dimensional  particles,  h,  w,  and  q  were  varied  to  independently  change 
the  dipolar  and  octopolar  moments  (from  symmetry  the  quadrupole  moments  are 
identically  zero),  and  mimic  the  interactions  of  the  magnetic  spheres  in  the  experi¬ 
ment.  This  model  reduces  to  a  two-charge  model  in  the  limit  of  w/D  — >  0.  ft  can 
also  be  reduced  to  a  two-point-dipole  model  in  the  limit  of  of  h/D  — >  0,  and  a  sin¬ 
gle,  centrally  located  point  dipole  model  if  h/D  ^  0  and  w/D  — >  0.  Similar  point 
charge  models  were  also  developed  to  describe  interactions  between  anisotropically 
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Fig.  2.3:  A  schematic  of  our  general  magnetic  particle  model.  The  model  consists  of  four 
point  charges  of  equal  magnitude  ±q,  separated  by  distances  h  and  w  (0  <  h, 
w  <  D),  such  that  the  charges  are  at  the  corners  of  a  rectangle. 

interacting  biological  macromolecules.  Those  will  be  described  in  Chapters  6  and  7. 


and  p(x)  is  the  charge  density  at  position  x.  In  general,  the  l  =  0  term  in  the 
expansion  is  named  the  ‘monopolc’  term,  the  l  =  1  term  in  the  expansion  is  named 
the  ‘dipole’  term,  the  /  =  2  term  in  the  expansion  is  named  the  ‘quadrupole’  term, 
l  =  3  term  in  the  expansion  is  named  the  ‘octopole’  term,  and  so  on  for  higher  order 
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terms.  In  general,  we  may  also  speak  of  ‘multipole  moments’  to  characterise  the 
various  terms  in  the  multipole  expansion.  In  this  thesis,  we  define  these  moments 
as  follows: 


•  Monopole  Moment  (scalar), 

N 

Q  =  ^2  Qi 

i=  1 


•  Dipole  Moments  (vector), 


N 

Pj  =  Qdi 

i=  1 

•  Quadrupole  Moments  (traceless  rank  2  tensor), 

N 

Qjk  'y  )  U i  ( '\}i P f  f . i  &jk) 

i=  1 


•  Octopole  Moments  (rank  3  tensor), 

N 

Ojkl  ^  )  Qijikili 
i— 1 


(2.3) 


(2.4) 


(2.5) 


(2.6) 


In  general,  these  moments  are  associated  with  their  respective  terms  in  the  expansion 
through  the  following  relation: 


$(x) 


i  Q  p  •  x 

47reo  r  r3 


Z  hJ 


(2.7) 
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3.  POLYMERIZATION  TRANSITIONS  IN  TWO-DIMENSIONAL 


SYSTEMS  OF  DIPOLAR  SPHERES 

(Simulations  by  Kevin  Van  Workum  and  Jack  F.  Douglas,  Polymers  Division,  Na¬ 
tional  Institute  of  Standards  and  Technology) 

3. 1  Overview 

In  this  chapter,  we  investigate  the  clustering  transition  in  our  two-dimensional  sys¬ 
tem  of  magnetic  spheres,  as  a  function  of  vibration  amplitude.  Following  previous 
work  on  vertically  vibrated  granular  matter  [44],  we  draw  an  analogy  between  vi¬ 
brational  excitation  in  our  experimental  system  and  thermal  excitation  in  a  Monte 
Carlo  simulation  of  a  similar  system.  We  show,  using  our  experimental  system  and 
Monte  Carlo  simulations  performed  by  collaborators  at  NIST ,  that  at  temperatures 
above  the  clustering  transition,  the  particles  are  largely  dissociated  from  each  other 
and  exhibit  traditional  gas  phase  characteristics.  At  temperatures  below  the  poly¬ 
merization  transition,  the  thermal  (kinetic)  energy  of  the  particles  is  insufficient  to 
overcome  the  attractive  interactions  and  the  particles  form  strongly  bound  clusters. 
The  anisotropic  nature  of  the  dipolar  interactions  induce  head-to-tail  alignment 
resulting  in  the  formation  of  linear,  ring,  or  branched  structures.  Between  these 


two  limits  lies  a  broad  transition  region  in  which  these  polymer-like  structures  con¬ 
tinuously  form  and  disappear  in  dynamic  equilibrium.  Our  results  show  that  the 
transition  temperature,  T$,  increases  and  the  width  of  the  transition  broadens  with 
increasing  concentration,  C  (measured  as  particle  area  fraction  in  our  2D  system), 
in  accord  with  the  analytic  theory  of  equilibrium  polymerization  [48]. 

It  is  important  to  have  well-defined  models  of  self-organization  that  contain 
the  minimal  physical  characteristics  of  those  systems  and  to  extensively  study  the 
essential  characteristics  of  these  systems.  We  need  to  determine  what  kind  of  ther¬ 
modynamic  transitions  are  involved  and  how  to  determine  the  relevant  interaction 
parameters.  The  dynamical  properties  of  these  systems  (formation  and  the  prop¬ 
erties  while  functioning)  are  also  a  crucial  concern  and  we  need  to  develop  simple 
models  that  are  amenable  to  large  scale  simulations  so  that  these  properties  can 
be  studied  with  existing  computational  resources.  Optimally,  the  physical  systems 
chosen  should  have  well-defined  physical  realizations  that  can  be  compared  to  both 
simulation  and  analytic  theories  of  self-organization. 

Van  Workum  and  Douglas  [49]  have  advocated  dipolar  fluids  as  a  prototypical 
model  of  self-organization  and  they  have  extensively  investigated  the  polymerization 
transition  which  occurs  in  this  model  by  Monte  Carlo  simulations  in  three  dimen¬ 
sions.  Stambaugh  and  Losert  [41,  42]  have  experimentally  and  theoretically  studied 
the  self-assembly  of  spheres  with  multipolar  interactions  and  have  found  sensitive 
dependance  of  the  self-assembly  pattern  on  higher  order  moments.  This  chapter 
considers  a  physical  realization  of  a  dipolar  fluid  constrained  to  be  nearly  two- 
dimensional,  and  a  corresponding  Monte  Carlo  simulation  of  this  two-dimensional 
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fluid.  This  comparison  not  only  provides  a  check  on  the  physical  model  developed 
involving  magnetic  particles  shaken  vertically  to  thermalize  the  granular  fluid,  but 
provides  a  valuable  check  on  the  thermodynamic  characteristics  provided  by  the 
Monte  Carlo  simulations  and  the  analytic  theory  of  equilibrium  polymerization. 
The  choice  of  the  dipolar  fluid  is  independently  motivated  by  the  fact  that  protein 
molecules  characteristically  have  large  electrostatic  dipole  moments  (e.g.  tubulin  « 
1400  Debye,  as  will  be  shown  in  Chapter  6)  and  dipolar  interactions  clearly  have 
relevance  to  the  self-organization  processes  occuring  in  proteins  [77]. 

3.2  Experimental  Details 

The  magnetic  particles  used  in  this  chapter  had  magnetic  cores  of  length,  l  =  1.42  ± 
0.01  cm.,  and  diameter,  d  =  0.94  ±  0.01  cm.  The  particles  ( long  particles)  were 
5.3  ±0.1  grams  in  mass,  the  maximum  strength  of  the  magnetic  field  on  the  surface 
of  a  particle  was  0.54  kG  ±16%,  and  the  particle  diameter,  D,  was  1.69  ±  0.01 
cm.  The  container  is  constructed  with  a  black  Delrin  bottom  and  N  particles  were 
placed  in  the  container,  to  give  a  system  concentration,  C  =  prior  to  vibration. 
The  system  was  imaged  from  above  at  a  rate  of  0.5  Hz  with  a  1280x1024  pixel  color 
camera,  a  was  varied  from  5 g  to  11  g. 

Earlier  work  has  shown  that  a  direct  analogy  can  be  made  between  the  excita¬ 
tion  of  macroscopic  particulate  matter  by  vertical  shaking  and  thermal  excitations 
of  microscopic  matter  [44] .  For  the  range  of  driving  amplitudes  studied,  this  analogy 
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can  be  quantified  by  the  expression 


T  oc  nj,  (3.1) 

where  T  is  the  effective  temperature  and  vp  is  the  peak  velocity  of  the  vertically 
vibrated  surface.  vp  can  be  easily  calculated  from  the  vibration  acceleration  and 
frequency,  and  hence,  we  define  our  experimental  temperature  as  T  =  vp. 


3.3  Simulation  Details 


Simulations  performed  by  Kevin  Van  Workum  and  Jack  F.  Douglas,  Polymers  Di¬ 
vision,  National  Institute  of  Standards  and  Technology. 

The  magnetic  particles  are  modeled  as  hard  spheres,  as  described  in  Chapter 
2,  but  with  w  =  0  such  that  they  are  made  up  of  two  oppositely  charged  monopoles. 
The  monopoles  are  separated  by  a  distance  h  =  0.69.D  and  are  situated  symmetri¬ 
cally  about  the  center  of  the  sphere.  The  potential  energy  between  monopole  i  on 
sphere  a  and  monopole  j  on  sphere  b  is  given  by 


u 


CLibj 


Qai  Qbj 


(3.2) 


where  rai&.  is  their  separation  distance  and  the  g’s  denote  the  charges  of  the  monopoles 
This  model  is  chosen  to  effectively  represent  the  real  interactions  of  the  magnetic 
particles  used  in  the  experiments. 

In  order  to  mimic  the  quasi- two-dimensional  nature  of  the  experimental  setup, 
the  centers  of  the  spheres  are  restricted  to  lie  in  a  plane  while  the  embedded  dipole 
vectors  are  allowed  to  rotate  in  3D.  All  the  spheres  also  interact  with  a  boundary 
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consisting  of  a  radially  symmetric  hard  wall  (a  cylinder)  similar  to  that  in  the 
experiments. 

Canonical  ensemble  Monte  Carlo  (MC)  simulations  are  used  to  generate  equi¬ 
librium  configurations  and  thermodynamic  averages.  In  addition  to  the  traditional 
translational  and  rotational  MC  moves,  two  advanced  MC  methods  are  imple¬ 
mented.  The  Aggregate  Volume  Bias  MC  method  [45]  is  used  to  enhance  the  dy¬ 
namic  formation  and  destruction  of  clusters.  This  method  provides  efficient  phase 
space  sampling  in  systems  that  form  strong  interparticle  associations  such  as  those 
present  in  this  study.  The  Parallel  Tempering  MC  method  [46,  47]  is  also  used 
in  these  simulations  to  provide  improved  statistical  sampling.  In  this  algorithm, 
several  independent  simulations  proceed  in  parallel,  each  at  a  different  temperature 
but  with  the  same  number  of  particles.  Periodically,  two  simulations  adjacent  to 
one  another  with  respect  to  their  temperatures  are  chosen  at  random  and  a  swap 
of  their  particle  configurations  is  attempted.  Attempted  swaps  are  accepted  with  a 
probability  given  by 

Pa^b  =  min[l,exp{At/A/3}],  (3.3) 

where  A U  =  Uf,  —  Ua  is  the  potential  energy  difference  of  the  systems,  and  A/3  = 
1  /(ksTb)  —  1  /(ksTa).  T,  for  the  simulations,  is  defined  in  units  of  the  dipolar  par¬ 
ticle  self-energy  divided  by  ks ■  At  low  temperatures,  the  particles  can  form  large, 
strongly  associated  clusters  that  may  become  trapped  in  local  potential  energy  min¬ 
ima  causing  the  exploration  of  relevant  phase-space  to  proceed  inefficiently.  At  high 
temperatures  the  particles  are  largely  free  from  such  constraints  and  the  simulation 
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samples  configuration  space  more  effectively.  By  ‘connecting’  the  simulations  at 
low  temperature  to  those  at  high  temperature  through  Parallel  Tempering,  one  im¬ 
proves  the  statistical  sampling  of  phase-space.  In  the  present  study,  24  simulations 
were  used  at  each  concentration  and  adjacent  simulations  differed  in  temperature 
by  AT  =  0.01. 


3.4  Results 

Based  on  previous  analytical  modeling  [48]  and  Monte  Carlo  simulations  [49]  of 
equilibrium  polymerization  of  the  Stockmayer  fluid  (dipolar  spheres  with  additional 
isotropic  Lennard- Jones  van  der  Waals  interactions),  we  anticipate  a  polymerization 
transition  to  occur  at  constant  concentration.  At  temperatures  above  the  polymer¬ 
ization  transition,  the  particles  are  largely  dissociated  from  each  other  and  exhibit 
traditional  gas  phase  characteristics.  At  temperatures  below  the  polymerization 
transition,  the  thermal  (kinetic)  energy  of  the  particles  is  insufficient  to  overcome 
the  attractive  interactions  and  the  particles  form  strongly  bound  clusters.  The 
anisotropic  nature  of  the  dipolar  interactions  induces  head-to-tail  alignment,  result¬ 
ing  in  the  formation  of  linear,  ring,  or  branched  structures.  Between  these  two 
limits  lies  a  broad  transition  region  in  which  these  polymer-like  structures  continu¬ 
ously  form  and  disappear  in  dynamic  equilibrium.  For  the  present,  we  loosely  define 
the  polymerization  transition  temperature,  T$,  as  the  temperature  at  which  the 
particles  first  begin  to  form  clusters.  A  more  rigorous  definition  of  based  on  the 
average  cluster  size  follows. 
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As  examples  of  the  observed  structures,  Fig.  3.1  show  the  experiment  (a-c) 
for  C  =  0.026  at  T/T$  >  1,  T/T$  ~  1  ,  and  T/T$  <  1,  respectively.  Fig.  3.1  (d-f) 
shows  similar  renderings  of  the  simulations  for  C  =  0.026  at  T/T$  >  1,  T /T$  «  1 
,  and  T/T<j>  <  1,  respectively.  For  T/T$  <  1,  a  linear  chain  of  head-to-tail  aligned 
particles  is  observed  and  for  T /T, j>  >  1  the  particles  are  completely  dissociated.  Near 
the  transition,  T/T$  w  1,  a  fraction  of  particles  are  loosely  associated  while  the 
remaining  are  dissociated.  This  is  onr  first  hint  that  the  observed  broad  transition 
involves  the  dynamic  formation  and  destruction  of  polymeric  chains  at  equilibrium 
with  larger  chains  forming  as  the  temperatures  decreases. 

We  first  characterize  the  transition  by  monitoring  the  average  cluster  size  as  a 
function  of  temperature.  Two  particles  are  considered  to  be  part  of  the  same  cluster 
when  their  interparticle  distance  is  less  than  ra  =  1.125a.  Based  on  previous  work 
[50],  the  observed  average  cluster  size  is  expected  to  be  relatively  independent  of 
small  changes  in  ra.  The  average  cluster  size  is  then  defined  by 

T  _  gjk  m 

Eh  N. \  ’ 

where  i  is  the  number  of  particles  in  a  given  cluster  and  is  the  number  of  chains 
of  length  i.  In  Fig.  3.2,  we  show  L  as  a  function  of  T  for  the  experiments  (a)  and  the 
simulations  (b)  for  a  range  of  concentrations.  We  find  a  family  of  curves  describing 
the  general  increase  of  L  upon  cooling.  At  higher  concentrations,  the  increase  in 
L{T)  occurs  faster  as  T  is  lowered,  reflecting  the  concentration  dependence  of  the 
polymerization  transition.  The  curves  in  Fig.  3.2  notably  have  a  similar  shape,  and 
it  is  natural  to  seek  a  reduced  variable  description. 
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Fig.  3.1:  A  photograph  of  an  experimental  system  where  the  concentration  C  =  0.026  at 
three  different  temperatures:  T/T$  >  1  (a),  T/T$  ~  1  (b),  and  T/T$  <  1  (c).  A 
rendering  of  the  simulation  where  C  =  0.026  and:  T /T$  >  1  (d),  T/T$  «  1  (e), 
and  T/T$  <  1  (f). 
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Fig.  3.2:  a)  A  plot  of  L  as  a  function  of  T  for  experimental  systems  where  C  =  0.020, 
0.052,  and  0.104.  b)  A  plot  L  as  a  function  of  T  for  similar  simulated  systems. 
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Previous  simulations  of  the  polymerization  transition  in  the  Stockmayer  fluid 
have  demonstrated  that  a  universal  reduced  variable  description  of  the  temperature 
dependence  of  L  can  be  obtained  by  expressing  temperature  relative  to  its  value 
of  the  polymerization  transition,  7$  defined  by  an  inflection  point  in  the  order 
parameter  for  the  polymerizatino  transition  <f>  (See  Eqn.  3.5  for  definition).  It  was 
also  found  that  L  at  this  temperature  is  nearly  universal,  taking  a  value  2.1  regardless 
of  particle  concentration  [49].  This  result  is  predicted  also  by  the  analytic  theory  of 
equilibrium  polymerization  [48].  We  checked  this  relation  for  the  simulated  analog 
of  onr  experimental  system  and  found  that  L  at  the  polymerization  transition  (7$) 
again  lies  in  a  narrow  range  for  the  concentrations  we  consider.  This  universality 
in  the  magnitude  of  L  at  7$  suggest  that  we  should  similarly  be  able  to  reduce  the 
scatter  in  the  measurements  in  Fig.  3.2  by  simply  normalizing  temperature  by  the 
temperature  where  L  —  2,  to  obtain  an  approximate  equation  of  state  description 
of  the  T  dependence  of  L.  This  procedure  is  motivated  by  the  comparitively  noisy 
nature  of  our  <h  data. 

We  see  from  Fig.  3.3(a)  that  this  procedure  indeed  reduces  the  scatter  con¬ 
siderably  and  the  reduction  is  quite  good  for  the  approach  to  the  polymerization 
transition  temperature  T  — ■»  7$  where  L  —  2.  Note  that  the  simulation  data  in 
Fig.  3.3(b)  exhibits  a  tendancy  to  saturate  to  a  finite  value,  which  is  simply  the 
number  of  particles  in  the  system.  This  feature  is  a  finite  size  effect  and  the  growth 
of  L  at  low  T  is  apparently  unbounded  in  the  thermodynamic  limit,  as  found  before 
for  the  three-dimensional  Stockmayer  fluid.  The  chains  in  the  experiment  also  vi¬ 
sually  exhibit  a  saturation  of  chain  length  to  a  size  corresponding  to  the  number  of 
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particles  in  the  system  at  low  T,  but  the  particle  tracking  algoithm  can  not  follow 
these  large  clusters  reliably,  giving  rise  to  the  residual  scatter  in  Fig.  3.3(a)  and 
the  absence  of  a  discernablc  plateau.  Nevertheless,  Fig.  3.3(a)  clearly  indicates  the 
sharp  rise  in  L  at  lower  T  and  the  significance  of  L(T$)  =  2  in  defining  a  reduced 
variable  description  of  these  observations.  We  also  note  that  our  determination  of 
<h(T)  below  is  consistent  with  L{T<$>)  ~  2  within  experimental  uncertainty  in  our 
experimental  system  and  we  could  equally  as  well  have  reduced  our  L  data  by  7$ 
values  determined  directly  fomr  the  inflection  point  of  <3>  Both  procedures  lead  to  a 
regettable  large  scatter  at  low  T,  where  the  uncertainties  from  our  particle  tracking 
estimates  of  L  and  are  large. 

It  is  also  natural  to  monitor  the  average  fraction  of  particles  that  are  members 
of  a  cluster,  i.e.  the  extent  of  polymerization,  $.  The  extent  of  polymerization  is 
defined  by 


$  = 


(3.5) 


where  Np  is  number  of  polymerized  (reacted)  particles  and  N  is  the  total  number 
of  particles  in  the  system.  therefore  defines  an  order  parameter  where  $  =  0 
corresponds  to  an  entirely  dissociated  system,  and  <3>  =  1  corresponds  to  a  state 
where  all  the  particles  are  in  contact  with  at  least  one  other  particle. 

Figure  3.4(a)  is  a  plot  of  $  versus  T  for  experimental  systems  where  C  = 
0.020,  0.052,  and  0.104.  A  transition  from  $  «  1  to  $  ~  0  is  seen  with  increasing 
T  for  C  =  0.020,  while  the  transitions  for  C  >  0.052  are  partially  captured  for  the 
temperature  range  studied.  While  it  is  likely  that  observations  of  the  C  >  0.052 
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Fig.  3.3:  a)  A  plot  of  L  as  a  function  of  T/T$  for  the  experimental  systems  where  C  = 
0.020,  0.052,  and  0.104.  b)  A  plot  of  L  as  a  function  of  T /T$  for  similar  simulated 
systems. 
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Fig.  3.4:  (a)A  plot  of  $  versus  T  for  experimental  systems  where  C  =  0.020,  0.052,  and 
0.104.  (b)  A  similar  plot  of  $  versus  T  for  similar  simulated  systems.  Note  that 
for  both  simulation  and  experiment,  T, j>  increases  and  the  transitions  broaden 
with  increasing  C. 
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systems  at  higher  T  would  show  transitions  similar  to  those  at  C  <  0.052,  the 
experimental  setup  limited  our  observations  to  T  <  0.3.  Figure  3.4(b)  is  a  plot  of  <f> 
versus  T  for  simulated  systems  for  C  =  0.020,  0.052,  and  0.104.  A  transition  from 
pa  1  to  <f>  ~  0  is  seen  with  increasing  T  for  all  C.  While  the  simulation  methods 
exclude  the  possibility  of  investigating  temporal  dynamics,  the  experimental  system 
allows  a  detailed  study  of  the  evolution  of  the  spacial  patterns  in  time. 

Figure  3.5  is  a  plot  of  $  versus  time  for  the  experimental  system  where  C  = 
0.052  and  T  «  T$.  varies  significantly  within  the  timespan  of  the  experiment, 
with  <f>  «  0  and  $  «  1  at  various  times,  though  T  remains  constant.  This  large 
variance  and  intermittency  is  characteristic  of  the  system  when  T  «  X$,  and  similar 
phenomena  are  observed  at  other  concentrations. 

Figure  3.6  is  the  distribution  of  <F’s  observed  in  Fig.  3.5,  normalized  to  the 
total  number  of  frames  in  the  experiment.  The  broad  range  of  the  distribution  is  a 
characteristic  of  the  system  at  T  T$. 

Figure  3.7  is  a  plot  of  the  variance  in  <f>  versus  T/T$  for  the  experimental 
system  where  C  =  0.052.  Note  the  significant  peak  at  T /T$  «  1.  This  peak  at 
T  «  T<j,  is  a  characteristic  behavior  of  the  system  at  the  transition,  and  it  is  also 
observed  at  other  concentrations. 

In  order  to  directly  compare  in  our  experiments  and  simulations,  we  rescale 
our  experimental  temperature  using  the  expression 

T  =  AT  +  B ,  (3.6) 


where  A  and  B  are  constants  to  be  determined.  This  method  of  rescaling  is  con- 


Fig.  3.5:  A  plot  of  <5  versus  time  for  the  experimental  system  where  C  =  0.052  and  T  ~  T$. 

Note  the  large  variation  in  <h,  with  $  ~  0  and  $  ~  1  at  various  times,  though  T 
remains  constant.  This  large  variance  and  intermittency  is  characteristic  of  the 
system  when  T  ~  T$. 
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Fig.  3.6:  A  plot  of  the  distribution  of  <f>’s  observed  in  Fig.  3.5,  normalized  to  the  total 
number  of  frames  in  the  experiment.  Note  that  the  distribution  is  nonzero  nearly 
over  the  entire  range  of  $. 
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Fig.  3. 7:  A  plot  of  the  variance  in  $  versus  T /T$  for  the  experimental  system  where  the 
concentration  C  =  0.052.  Note  the  significant  peak  at  T/T$  ~  1. 


sistent  with  the  thermal  relation  for  vertically  vibrated  systems  described  in  earlier 
work  [44].  A  and  B  are  determined  from  a  linear  best  fit  to  a  plot  of  the  7$  values 
from  the  experiment  (at  C  =  0.020,  0.026,  0.039,  0.052,  0.078,  and  0.104)  vs.  the 
7$  values  from  the  simulation  at  similar  C .  For  our  experiments,  A  =  13.798  and 
B  =  0.15144.  Figure  3.8  shows  the  polymerization  transition  temperatures,  7$  and 
7$  (  7$  are  the  rescaled  experimental  7$  values),  as  a  function  C  for  the  experi¬ 
ment  (filled  circles)  and  the  simulation  (open  circles).  Also,  in  Fig.  3.8  is  a  fit  of  the 
Dainton-Ivin  equation  [51]  to  the  simulation  data, 


n  = 


A  sp  +  kslnC ’ 


where  A hv  =  —2.13  ±  0.19,  and  A Sp/ks  =  —8.5  ±  1.0  (in  units  of  particle  self 


energy).  A  similar  fit  to  the  experimental  data  gives  A hp  =  —2.20  ±  0.34,  and 
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Fig.  3.8:  (open  circles)  and  (filled  circles)  versus  the  concentration  C  for  the  systems 

studied.  Also  plotted  is  a  fit  of  to  the  Dainton-Ivin  equation. 

A sp/kB  =  —8.9  ±  1.8.  Note  that  7$  and  increase  with  increasing  C ,  in  a  similar 
fashion  for  both  the  experiment  and  the  simulation. 

3.5  Conclusion 

We  have  investigated  the  self-organization  of  dipolar  spheres  into  chains  as  a  fun¬ 
damental  model  of  the  self-assembly  of  particles  having  anisotropic  interparticle 
interactions.  Our  studies  of  a  two-dimensional  system  of  confined  dipolar  spheres 
has  shown  that  there  exists  a  polymerization  transition  of  the  spheres  at  finite  tem¬ 
peratures.  Using  a  simple  model  system  of  vertically  vibrated  magnetic  spheres, 
together  with  Monte  Carlo  simulations  of  a  similar  system,  carried  out  by  collabo¬ 
rators  at  NIST,  we  showed  that  the  transition  temperature  and  the  broadness  of  the 
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transition  is  concentration  dependent.  Our  results  show  that  the  transition  tem¬ 
perature,  T$,  increases  and  the  width  of  the  transition  broadens  with  increasing 
concentration,  C,  in  accord  with  the  analytic  theory  of  equilibrium  polymerization 
of  Dudowicz  et  al.  [48] . 


3.6  Outlook 

Possible  future  work  includes  investigating  the  role  of  the  system  size  in  the  tran¬ 
sitions.  This  investigation  could  be  done  experimentally  by  building  containers  of 
different  diameters  and  repeating  the  experiments  at  similar  concentrations  as  those 
investigated  above.  Such  studies  would  help  determine  what  role  the  system  size 
has  on  our  results,  and  allow  for  better  statistics  in  the  case  of  larger  systems. 

It  would  also  be  interesting  to  investigate  the  role  of  higher  order  interactions 
on  these  transitions.  Quadrupolar  particles  could  be  built,  as  discussed  in  Chapter 
2,  or  the  role  of  the  octopolc  could  be  investigated  by  using  cylindrical  magnetic 
cores  of  differing  lengths.  These  studied  would  give  insight  into  the  affects  of  higher 
order  interactions  on  polymerization  transitions,  and  may  have  specific  implications 
for  biological  macromolecular  self-assembly  where  higher  order  interactions  are  es¬ 
sential. 

Another  possible  extension  of  this  work  is  to  build  sturdier  magnetic  particles 
with  weaker  interactions  so  that  higher  relative  temperatures  can  be  investigated 
without  destroying  the  particles. 
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4.  PATTERN  FORMATION  IN  A  MONOLAYER  OF  MAGNETIC 


SPHERES 

4.1  Overview 

In  this  chapter,  we  report  observations  of  pattern  formation  in  our  experimental  sys¬ 
tem  which  depends  on  the  shape  of  the  magnetic  core,  i.e.  the  magnetic  multipole 
moment  distribution.  While  systems  of  dipolar  spheres  have  been  investigated  rela¬ 
tively  extensively  [6,  7,  8,  9,  10,  11,  12,  26,  29,  30,  31,  34,  43],  the  roles  of  higher  or¬ 
der  moments  in  self-assembly  are  relatively  unknown.  Our  results  show  that  spheres 
containing  short  cylindrical  magnets  form  HOP  (hexagonal  close-packed)  maerovor- 
tices(Fig.  4.1(a)),  while  spheres  containing  long  cylindrical  magnets  form  concentric 
rings.  This  implies  that  the  pattern  formation  depends  not  only  on  the  dipolar  mo¬ 
ment,  the  temperature  of  the  system,  and  the  concentration  of  the  spheres,  but  also 
on  higher  order  multipole  moments.  In  particular,  because  of  the  cylindrical  sym¬ 
metry  of  the  magnetic  cores,  our  work  illustrates  the  effect  of  the  octopole  moment 
on  the  self-assembly. 

Monte  Carlo  simulations  of  the  thermal  equilibrium  of  large  2D  systems  of 
hard  spheres  with  centrally  embedded  point  dipoles  have  shown  that  ordering  in 
the  absence  of  external  forces  depends  on  the  dipolar  strength,  the  temperature  of 


Fig.  4.1:  (a)  A  schematic  illustration  of  the  ‘HCP  macrovortex’  pattern.  The  dipoles  are 
oriented  in  a  macroscopic  circulating  vortical  pattern  about  the  center,  (b)  A 
schematic  illustration  of  the  ‘square-packed  microvortex’  pattern.  The  dipole 
axes  are  tilted  by  about  45°  with  respect  to  the  axes  of  the  square  pattern, 
and  the  dipole  orientations  are  such  that  they  define  alternating  directions  of 
circulation  (indicated  by  the  rotational  arrows  in  the  centers  of  the  squares)  on 
a  checker-board  pattern  of  squares  (after  Ref.  [52]). 

the  system,  and  the  number  density  of  the  spheres  [11,  29,  30,  31].  At  interme¬ 
diate  temperature,  high  density  systems  of  dipolar  spheres  form  hexagonal-close- 
packed  (HCP)  macrovortices  (Fig.  4.1(a)),  while  intermediate  density  systems  form 
networks  of  long  branching  chains  [29].  At  high  temperatures,  systems  of  dipolar 
spheres  disassociate  and  become  gaseous  [26,  29]. 

A  calculation  done  by  Belobrov  et  al.  [52]  showed  that  the  zero-temperature 
ground  state  of  a  finite  (<  500  particle)  two-dimensional  HCP  lattice  of  dipoles  is  a 
macrovortex  (Fig.  4.1(a)),  while  the  zero-temperature  ground  state  of  a  finite  square- 
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packed  lattice  of  dipoles  is  a  microvortex  (Fig.  4.1(b)).  The  microvortex  state  as  a 
minimum  of  the  interaction  Hamiltonian,  was  also  remarked  upon  by  Rozenbaum  et 
al.  [53].  A  similar  macrovortex  state  seen  at  the  atomic  scale  in  500  nm  diameter  thin 
magnetic  films  is  also  of  significant  current  technological  interest  as  a  novel  magnetic 
RAM  material  [16,  17,  18,  19,  20,  21,  22,  23].  The  possibility  of  information  storage 
in  such  systems  occurs  due  to  an  ordered  out  of  plane  component  of  the  atomic 
spins  near  the  center  of  the  macrovortex.  While  macrovortecies  have  been  observed 
in  magnetic  thin  films  of  thickness  less  than  30  nm,  the  out  of  plane  ordering  has 
been  shown  to  vanish  in  the  case  of  magnetic  monolayers.  The  central  out  of  plane 
dipole  orientation  is  independent  of  the  ordering  of  the  rest  of  the  macrovortex,  and 
is  easily  flipped  electronically,  allowing  for  the  possibility  of  two  bits  of  memory 
storage  for  each  micron-diameter  disk.  The  macrovortex  state  was  also  observed 
in  simulations  of  a  dense  system  of  dipolar  spheres  in  thermal  equilibrium  [29], 
and  in  an  experimental  system  of  settled  magnetic  spheres  (a  system  of  magnetic 
spheres  that  are  initially  magnetically  excited  with  a  solenoid  and  then  allowed 
to  settle  at  the  bottom  of  a  container  under  zero  external  magnetic  field)  [15]. 
However,  to  our  knowledge,  a  macrovortex  pattern  has  never  been  observed  in  a  finite 
temperature  experimental  system  of  dipolar  spheres.  We  provide  the  first  direct 
experimental  demonstration  of  the  macrovortex  pattern  in  a  thermalized  system  of 
dipolar  spheres.  We  also  provide  a  direct  experimental  visualization  of  a  microvortex 
pattern  in  a  zero  temperature  system  of  dipolar  spheres. 
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4.2  Experimental  Details 


In  the  following  experiments,  we  investigate  particles  encasing  magnetic  cores  of 
two  types  which  we  call  ‘long’  and  ‘short’  magnets.  The  long  magnets  we  used 
had  length,  l  =  1.42  ±  0.01  cm.,  and  diameter,  d  =  0.94  ±  0.01  cm.  The  short 
magnets  we  used  had  l  =  0.64  ±  0.01  cm.  and  d  =  0.95  ±  0.01  cm.  Once  the  long 
magnets  had  been  encapsulated  in  the  spherical  shells,  the  particles  ( long  particles) 
were  5.3  ±  0.1  grams  in  mass,  the  maximum  strength  of  the  magnetic  field  on  the 
surface  of  a  long  particle  was  0.54  kG  ±16%.  The  particles  containing  the  short 
magnets  ( short  particles)  were  4.3  ±  0.1  grams  in  mass,  the  maximum  strength  of 
the  magnetic  held  on  the  surface  of  a  short  particle  was  0.12  kG  ±24%.  The  short 
magnets  were  secured  inside  the  plastic  shells  by  filling  the  extra  space  in  each  half 
shell  with  RTV  resin.  All  particles  were  constructed  to  have  a  light  colored  north 
hemisphere  and  a  dark  colored  south  hemisphere  so  that  dipole  orientations  could  be 
discerned.  (Previous  experiments  on  thermally  or  externally  excited  systems  did  not 
sense  the  dipole  orientations  of  individual  particles.)  The  container  was  constructed 
with  a  polyvinyl  chloride  bottom  and  N  particles  were  placed  in  the  container  prior 
to  vibration,  to  give  a  system  concentration,  C,  as  defined  in  Chapter  3.  Typically, 
0.30  <  C  <  0.65.  The  system  was  imaged  at  a  rate  of  2  Hz  for  500  frames  with 
a  1024x1024  pixel  camera.  In  these  experiments,  typical  accelerations  for  the  long 
particles  were  4  —  8 g,  while  short  particles  were  usually  driven  at  about  2  —  5 g. 
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4.3  Results 


Figure  4.2(a)  shows  the  ‘network’  pattern  of  short  particles  where  C  =  0.309,  after 
being  shaken  at  4.3 g  for  250  seconds.  In  this  case  (Fig.  4.2(a)),  the  initial  state  was 
created  by  pouring  the  particles  randomly  into  the  container.  The  observed  chain- 
like  structure  is  similar  to  the  pattern  observed  by  Blair  et  al.  [26]  in  intermediate 
density  systems  of  vertically  vibrated  magnetic  spheres,  and  by  Weis  [29]  in  simu¬ 
lations  of  intermediate  density  systems  of  dipolar  spheres.  Figure  4.2(b)  shows  the 
configuration  of  short  particles  where  C  =  0.618  after  being  shaken  at  4.3ry  for  250 
seconds.  The  initial  state  of  the  particles  was  random.  The  final  state  of  the  particles 
after  being  driven  is  mostly  an  HCP  macrovortex  state  similar  to  the  ground  state 
of  a  dipolar  hexagonal  lattice  predicted  by  Belobrov  et  al.  [52]  (Fig.  4.1(a))  and  seen 
in  simulations  by  Weis  [29].  We  find  that  this  HCP-based  macrovortex  is  a  stable 
state  for  the  short  particles.  When  other  patterns  such  as  square-packed  microvor¬ 
tices  (Fig.  4.1(b),  4.2(c))  or  concentric  rings  are  set  up  as  the  initial  condition,  after 
driving,  the  system  always  becomes  an  HCP-based  macrovortex  (Fig.  4.2(d)).  If  the 
system  is  set  up  as  a  macrovortex,  it  remains  in  a  macrovortex  state.  Also,  when 
several  single  chains  of  the  particles  are  brought  together  in  an  antiferromagnetic 
square-packed  state  of  190  particles  (each  chain  has  the  opposite  dipolar  orientation 
from  that  of  the  two  neighboring  chains),  they  immediately,  without  driving,  reori¬ 
ent  themselves  into  the  microvortex  state  (Fig.  4.1(b),  4.2(c)).  If  the  particles  are 
initially  set  up  in  a  ferromagnetic  HCP  pattern,  under  small  (2.4g)  agitation,  the 
macrovortex  state  is  assumed. 
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Fig.  4.2:  (a)  The  network  pattern  of  an  intermediate  density  (C  =  0.309)  system  of  short 
particles  after  being  shaken  at  4.3 g  for  250  seconds,  (b)  The  final  macrovortex 
configuration  of  a  dense  {C  =  0.618)  system  of  short  particles  after  being  shaken 
at  4.3(7.  (c)  The  nricrovortex  pattern  made  by  the  square-packed  set  up,  un¬ 
driven,  dense  (C  =  0.618)  system  of  short  particles,  (d)  A  ‘hollow  macrovortex’ 
configuration  of  the  short  particles  after  being  set  up  in  an  aligned,  concentric 
rings  pattern  and  then  shaken  at  4.3 g. 
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Fig.  4.3:  (a)  The  final  concentric  ring  state  of  a  dense,  C  =  0.618,  system  of  initially 
random  long  particles  after  being  shaken  at  7.8 g  for  250  seconds,  (b)  The  final 
state  of  the  system  after  being  shaken  in  a  square  box  at  7.8(7  f°r  250  seconds. 
The  pattern  is  one  of  squared-off  concentric  rings,  similar  to  the  final  state  of 
the  particles  in  the  cylindrical  container  of  equal  height  and  area. 

Figure  4.3(a)  is  the  final  concentric  ring  state  of  a  C  —  0.618  system  of  long 
particles  after  being  shaken  at  7.8 g.  The  initial  state  of  the  particles  was  random. 
The  existence  of  the  concentric  ring  state  is  independent  of  the  initial  conditions  of 
the  system.  If  the  system  is  set  np  in  any  other  state,  it  rearranges  into  a  concentric 
ring  state  with  sufficient  driving.  Systems  of  all  lower  densities  also  prefer  ring  and 
line  states  at  intermediate  and  low  accelerations  (<  8 g). 

Figure  4.4(a)  is  the  initial  random  state  of  a  dense  (190  particle)  system  of 
long  particles.  After  100  seconds  of  shaking  at  7.8 c/,  the  system  evolves  into  the 
spiral  pattern  shown  in  Fig.  4.4  (b)  (bine  and  red  lines  are  superimposed  upon  the 
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Fig.  4.4:  (a)  The  initial  random  state  of  a  dense  (190  particle)  system  of  long  particles. 

(b)  An  intermediate  state  of  the  same  system  of  long  particles  after  shaking  at 
7.8 g  for  100  seconds.  Note  the  spiral  ordering  that  occurs  (highlighted  by  blue 
and  red  drawn  lines),  (c)  The  final  state  of  the  system  after  250  seconds,  a  steady 
state  pattern  of  concentric  rings,  similar  to  that  shown  in  Fig.  4.3. 

pattern  to  highlight  the  spiral  ordering).  Figure  4.4(c)  is  the  final  state  of  the  system 
after  250  seconds,  a  steady  state  pattern  of  concentric  rings,  similar  to  that  shown 
in  Fig.  4.3. 

Figure  4.5  (a)  is  the  initial  state  of  an  intermediate  density  (C  =  0.325)  sys¬ 
tem  of  long  particles.  The  system  quickly  evolves  to  form  chains  of  variable  length 
(highlighted  by  red  drawn  lines)  after  30  seconds  of  vibration  at  6.4r/,  as  shown  in 
Fig.  4.5  (b).  Figure  4.5  (c)  is  the  final  state  of  the  system  after  100  seconds,  a  steady 
state  pattern  of  concentric  rings,  similar  to  that  shown  in  Fig.  4.3  and  Fig.  4.4  (c). 
Note  the  contrast  in  between  the  ring-like  ordering  of  the  long  particles  shown  in 
Fig.  4.5  (c),  and  the  network  pattern  formed  by  the  system  of  short  particles  at 
similar  density  shown  in  Fig.  4.2  (a).  The  ring  state  is  a  robust  steady  state  phe- 
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Fig.  4.5:  (a)  The  initial  state  of  an  intermediate  density  (C  =  0.325)  system  of  long 
particles,  (b)  An  intermediate  state  of  the  same  system  of  long  particles  after 
shaking  at  6.4 g  for  30  seconds.  Note  the  appearance  of  long  chains,  (highlighted 
by  red  drawn  lines),  (c)  The  final  state  of  the  system  after  100  seconds,  a 
steady  state  pattern  of  concentric  rings,  similar  to  that  shown  in  Fig.  4.3  and 
Fig.  4.4.  Note  that  the  ring  state  is  a  robust  steady  state  phenomenon  for  the 
long  particles,  reproducibly  appearing  at  several  different  particle  densities  and 
in  different  shaped  containers. 

nomenon  for  the  long  particles,  reproducibly  appearing  at  several  different  particle 
densities  and  in  different  shaped  containers,  as  illustrated  in  Figs.  4.3,  4.4,  and  4.5. 

We  can  track  the  positions  of  all  particles  during  our  experiments,  using  dig¬ 
ital  image  algorithms  developed  in  the  IDL  programming  environment.  To  gauge 
quantitatively  how  much  HCP  clustering  occurred  in  each  experiment,  we  designed 
a  program  that  finds  the  number  of  particles  with  6  touching  neighbors  in  each  im¬ 
age.  The  program  identifies  particles  as  ‘touching’  if  the  centers  of  the  particles  are 
within  9/8  particle  diameters  of  each  other.  We  found  that  for  a  system  of  190  short 
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particles  (C  =  0.618)  initially  in  a  random  state  and  then  shaken  at  4.3 g,  the  mean 
number  of  particles  with  6  touching  neighbors,  1V6  ,  averaged  over  500  frames,  was 
43.3  ±  9.0.  For  190  initially  random  long  particles  shaken  at  7.8g,  Nq  —  0.14  ±  .09. 
In  these  two  experiments,  where  190  particles  began  in  a  random  state,  the  N6  of 
the  short  particles  was  more  than  300  times  the  Nq  of  the  long  particles.  Using  this 
measure,  in  all  experiments  with  similar  initial  conditions,  the  short  particles  were 
quantitatively  more  clustered  than  the  long  particles. 

Boundary  conditions  do  not  appear  to  affect  the  qualitative  features  of  the  final 
states  for  either  type  of  particle.  We  performed  similar  experiments  in  a  square 
box  of  height  and  area  equal  to  those  of  the  cylindrical  container,  and  the  final 
states  were  always  qualitatively  similar  to  the  final  state  in  the  cylindrical  container. 
For  example,  see  Fig.  4.3(b),  which  shows  the  state  of  a  C  —  0.618  system  of 
long  particles  after  shaking  in  the  square  box  under  similar  initial  conditions  as  for 
Fig.  4.3(a). 


4.4  Energetic  Considerations 

To  gain  insight  into  the  reasons  why  the  long  particles  form  different  patterns  than 
the  short  particles,  we  numerically  examined  the  theoretical  model  described  in 
Chapter  2  with  w  — >  0.  This  model  is  similar  to  that  used  in  simulations  by 
Schneider  et  al.  [27]  to  model  lipid  headgroups  (which  are  also  dipolar),  and  it  gives 
qualitatively  similar  results  to  the  two-parallel-dipole  model  used  by  Camp  et  al.  [54], 
Similar  three  point  charge  models  are  used  to  match  charge,  dipole  and  quadrupole 
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moments  of  polar  molecular  systems  such  as  hydrogen  fluoride  for  simulations  [24], 
Our  model  also  reduces  to  a  central  point  dipole  model  as  h/D  — >  0. 

We  use  the  usual  Coulomb  energy  formula, 


_ Ho  mj 

47t  Rij 


(4.1) 


(similar  to  Eq.  3.2  in  Chapter  3)  and  sum  over  all  interparticle  pairs  of  charges  to 
obtain  the  energy  of  a  pattern  of  particles  (excluding  the  energy  of  pairs  of  charges 
within  the  same  particle).  When  we  take  the  limit  h/D  — >  0,  our  calculations  are 
consistent  with  the  energy  calculations  of  rings  and  lines  of  spheres  containing  single 
point  dipoles  done  by  Clarke  et  al.  [25]. 

It  is  interesting  to  note  that,  assuming  the  dipole  moments  are  aligned  with 
the  z  axis,  while  -P0.75/-P0.001  =  1  by  design  (where  Ph/d  is  the  z  component  of 
the  dipole  moment  of  a  model  particle  with  a  given  value  of  h/D ,  the  ratio  of  the 
zzz  components  of  the  octopole  moments,  O0.75/O0.001  =  562500.  Because  the  zzz 
component  of  the  octopole  moment  is  the  highest  order  differing  term,  it  is  likely 
the  most  important  differing  moment  for  models  of  different  h/D  for  small  h/D. 

In  Fig.  4.6,  we  have  plotted  the  energies  of  two  touching  particles  as  the  angle 
between  their  axes  is  symmetrically  bent  for  several  values  of  h/D.  Note  that  the 
slope  of  the  curves  increases  as  h/D  increases.  This  indicates  that  chains  of  particles 
with  large  h/D  are  stiffer  than  those  with  small  h/D.  Also  note  that  at  0  =  0,  the 
particles  are  in  the  head-to-tail  configuration  and  at  0  =  7t/2,  the  particles  are  in 
the  antialigned  side-to-side  configuration.  For  large  h/D ,  the  head-to-tail  (0  =  0) 
configuration  is  much  more  favorable  than  the  antialigned  side-to-side  (0  =  7t/2) 
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configuration  (e.g.  for  h/D  =  0.75,  the  head-to-tail  energy  is  more  than  6  times  the 
antialigned  side-to-side  energy).  For  h/D  =  0.001,  the  head-to-tail  energy  is  only 
twice  the  antialigned  side-to-side  energy.  This  fact  is  elucidated  by  plotting  the 
ratio  of  the  energy  of  two  particles  in  head-to-tail  alignment,  EHT  to  the  energy  of 
two  particles  in  an  antiparallel  side-to-side  alignment,  Ess,  versus  h/D ,  as  is  shown 
inf  Figure  4.7.  This  indicates  that,  as  the  distance  between  the  charges  increases, 
the  relative  importance  of  the  side-to-side  interaction  decreases.  Particles  with  large 
h/D  (e.g.,  our  long  particles)  have  a  strong  preference  to  align  head-to-tail,  while 
particles  with  small  h/D  (e.g.,  our  short  particles)  do  not.  This  is  in  contrast  to 
recent  models  of  dipolar  ’spherocylinders’  (cylinders  with  hemispherical  ends)  which 
encapsulate  a  single  centrally  located  point  dipole  [56,  57].  These  spherocylindrical 
particles  increase  the  relative  importance  of  side-to-side  interactions. 

Figure  4.8  is  a  plot  of  the  energy  gained  by  bringing  together  ten  10-particle 
chains  in  various  patterns  normalized  to  the  energy  of  ten  infinitely  separated  10- 
particle  chains  versus  h/D.  At  low  h/D  (<  0.19),  the  lowest  energy  state  is  the 
HCP  macrovortex  (as  predicted  by  Belobrov  et  al.  [52]  for  h/D  =  0  particles).  At 
h/D  >  0.19,  the  particles  energetically  prefer  a  ring  state.  This  is  consistent  with 
our  observations  of  high  density  systems;  particles  with  short  cylindrical  magnetic 
cores  form  a  macrovortex  (Fig.  4.2(b)),  while  particles  with  long  magnetic  cores 
form  rings  (Fig.  4.3(a)).  Also,  the  orientations  of  the  long  particles  in  each  ring 
appear  to  be  independent  of  the  orientations  of  the  particles  in  neighboring  rings. 
The  square-packed  microvortex  state  is  more  favorable  than  the  square-packed  an¬ 
tiferromagnetic  arrangement  only  at  very  small  h/D  (<  0.1).  This  is  also  consistent 


45 


Energy 


Fig.  4.6:  The  energies  in  units  of  of  two  touching  particles  as  a  function  of  the  angle 

0  between  their  axes  for  several  values  of  h/ D. 
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Fig.  4.7:  The  ratio  of  the  energy  of  two  particles  in  head-to-tail  alignment,  Eht  to  the 
energy  of  two  particles  in  an  antiparallel  side-to-side  alignment,  E$s,  versus  h/ D. 
Note  that  particles  with  large  h/ D  strongly  prefer  a  head-to-tail  alignment  (with 
Ejjt /Ess  ~  6  at  h/D  =  0.7),  while  particles  with  smaller  h/D  have  relatively 
little  preference. 
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with  our  observations.  When  we  attempt  to  set  np  the  short  particles  in  a  square- 
packed  antiferromagnetic  pattern,  they  immediately,  without  driving,  reorient  into 
a  microvortex  pattern.  The  long  particles  never  form  a  microvortex  pattern  when 
set  up  in  a  square-packed  antiferromagnetic  state,  even  when  driven.  The  HCP 
ferromagnetic  state  is  energetically  unfavorable  compared  to  separated  lines  at  all 
h/D.  This  is  in  contrast  to  the  result  for  an  infinite  HCP  lattice,  in  which  case 
the  ground  state  is  ferromagnetic  [58] .  It  is  also  interesting  to  note  that  the  ground 
state  of  an  HCP  lattice  of  magnetic  dipoles  can  change  under  an  external  tangential 
magnetic  field.  Bologa  et  al.  [15]  observed  a  transition  in  an  unexcited  system  of 
magnetic  spheres  from  an  HCP  macrovortex  state  to  an  HCP  ferromagnetic  state 
under  an  increasing  external  tangential  magnetic  field. 

We  conclude  that  particles  with  a  large  h/D ,  such  as  our  long  particles,  are 
more  likely  to  form  long,  stiff,  weakly  interacting  chains  and  rings.  Particles  with 
small  h/D  are  more  likely  to  aggregate  into  clustered  patterns.  Energetically,  it 
has  been  shown  that  the  ground  state  of  a  hexagonal  lattice  of  point  dipoles  is 
a  macrovortex  [52],  It  has  also  been  shown  that  for  systems  of  between  4  and  13 
(point)  dipolar  particles,  the  global  ground  state  is  a  single  ring  [25].  For  a  system  of 
14  or  more  dipolar  spheres,  a  ‘double  ring  state’  has  lower  energy,  and  it  is  expected 
that  the  lowest  energy  configuration  changes  once  again  to  triple,  quadruple,  and 
higher  order  ring  states  as  the  number  of  particles  increases.  We  conjecture  that,  as 
the  number  of  particles  increases,  the  global  ground  state  may  begin  to  look  like  a 
‘hollow  macrovortex’,  similar  to  that  in  Fig.  4.2(d). 

Figure  4.9  is  a  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely 
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h/D 


Fig.  4.8:  The  energy  of  ten  10-particle  chains  brought  together  in  various  patterns  nor¬ 
malized  to  the  energy  of  ten  infinitely  separated  10-particle  lines  versus  h/D. 
The  energy  of  two  infinitely  separated  10-particle  lines  is  normalized  to  0  for  all 
h/D.  Note  that  the  HCP  nracrovortex  is  the  lowest  energy  conformation  only  for 
h/D  <  0.19;  for  h/D  >  0.19,  a  single  100-particle  ring  state  is  of  lower  energy. 
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separated  chains  of  length  y/N)  of  several  different  patterns  of  dipolar  (h  =  0.001) 
particles  versus  y/N.  All  patterns  except  the  ring  pattern  are  rectangular  lattices 
of  dimensions  N  particles  by  N  particles.  Note  that  the  ring  state  is  energetically 
preferred  for  <  9,  while  for  y/N  >  9  the  HCP  macrovortex  is  preferred,  ft  is 
also  interesting  to  note  that  while  the  normalized  energy  of  the  HCP  macrovortex 
steadily  increases  for  y/N  >  20,  the  normalized  energy  of  the  HCP  ferromagnetic 
pattern  decreases.  By  taking  a  linear  best  fit  to  the  final  four  data  points  of  both  the 
HCP  patterns,  and  extrapolating  to  higher  y/N,  we  predict  the  ferromagnetic  state 
to  be  of  lower  energy  at  y/N  ps  105.  ft  is  well-known  that  there  is  a  transition  from 
macrovortex  to  ferromagnetic  ordering  in  an  HCP  dipolar  lattice  for  large  N  [52,  58], 
and  our  calculations  suggest  the  transition  would  occur  at  N  ps  11000. 

Figure  4.10is  a  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely 
separated  chains  of  length  v^A)  of  several  different  patterns  of  h  =  0.69  model 
particles  versus  y/N.  All  patterns  except  the  ring  pattern  are  rectangular  lattices 
of  dimensions  N  particles  by  N  particles.  Note  that  the  ring  state  is  energetically 
preferred  for  all  calculated  values  of  N  >2,  while  the  microvortex  and  macrovortex 
patterns  are  highly  unfavorable,  in  contrast  to  the  favored  patterns  of  the  h  =  0.001 
dipolar  particles  as  shown  in  Fig.  4.9. 

4.5  Conclusions 

Our  experimental  and  theoretical  results  show  that  the  shape  of  the  magnetic  core  in 
particles  can  be  a  determining  factor  in  pattern  formation,  even  in  simple  systems. 
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Normalized  Energy 


Fig.  4.9:  A  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely  separated  chains 
of  length  y/N )  of  several  different  patterns  of  dipolar  ( h  =  0.001)  particles  versus 
y/~N.  All  patterns  except  the  ring  pattern  are  rectangular  lattices  of  dimensions 
N  particles  by  N  particles.  Note  that  the  ring  state  is  energetically  preferred  for 
VN  <  9,  while  for  y/~N  >  9  the  HCP  nracrovortex  is  preferred. 
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Fig.  4.10:  A  plot  of  the  energies  (normalized  to  the  energy  of  \fN  infinitely  separated 
chains  of  length  \/N)  of  several  different  patterns  of  h  =  0.69  model  particles 
versus  y/N.  All  patterns  except  the  ring  pattern  are  rectangular  lattices  of 
dimensions  N  particles  by  N  particles.  Note  that  the  ring  state  is  energeti¬ 
cally  preferred  for  all  calculated  values  of  N  >  2,  while  the  nricrovortex  and 
macrovortex  patterns  are  highly  unfavorable,  in  contrast  to  the  favored  patterns 
of  the  h  =  0.001  dipolar  particles  as  shown  in  Fig.  4.9. 
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As  h/D  increases,  the  energetically  favored  pattern  changes  qualitatively.  Indeed, 
as  h/D  — >  1  in  our  model,  the  head-to-tail  arrangement  of  the  dipoles  becomes 
infinitely  more  favorable  than  any  other  configuration.  Our  unique  particles  allow 
us  to  provide  the  first  direct  experimental  evidence  of  the  square-packed  microvortex 
pattern  in  a  zero  temperature  system  of  dipolar  spheres,  and  the  first  observations 
of  an  HOP  macrovortex  state  in  an  excited  system  of  dipolar  spheres.  Finally,  we 
note  that  the  observed  patterns  can  be  understood  as  either  minimum  magnetostatic 
energy  equilibria  (in  the  macrovortex  and  microvortex  cases)  or  (in  the  case  of  low 
density  network  patterns,  Fig.  4.2(a))  as  thermal  equilibrium  states. 

4.6  Outlook 

Possible  future  work  includes  studies  of  larger  systems  of  particles  with  which  ori¬ 
entations  can  be  discerned.  This  could  be  achieved  by  building  larger,  lightweight 
containers  that  could  be  vibrated  on  a  shaker  similar  to  the  one  used  here.  Inves¬ 
tigations  of  larger  systems  would  allow  for  better  statistics,  but  more  importantly, 
it  is  possible  that  the  qualitative  features  of  the  observed  patterns  could  change  in 
large  systems.  Specifically,  for  dipolar  particles  (without  higher  order  interactions), 
it  has  been  shown  that  the  favored  pattern  of  an  HCP  lattice  is  ferromagnetic  for 
large  systems,  and  vortex  for  small  systems  [52,  58]. 

It  would  also  be  interesting  to  study  similar  systems  in  3D.  One  possible 
method  of  investigating  3D  systems  experimentally  would  be  to  build  colloidal  mag¬ 
netic  particles  with  variable  held  shapes,  and  density  match  them  to  a  surrounding 
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fluid.  Thermal  energy  could  then  provide  the  excitation  force,  rather  than  shaking. 
It  may  be  simpler  to  build  larger  density-matched  magnetic  particles,  and  then  ex¬ 
cite  them  with  an  AC  magnetic  field  imparted  by  a  soleniod,  similar  to  the  method 
used  by  Bologa  et  al.  [15]. 
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5.  SEGREGATION  IN  A  MONOLAYER  OF  MAGNETIC 


SPHERES 

5.1  Overview 

In  this  chapter,  we  show  using  our  experimental  system,  that  in  dilute  systems, 
the  long  particles  can  form  a  loose  unbranched  network  pattern.  If  the  multipole 
moment  distribution  of  a  small  fraction  of  the  magnetic  particles  is  altered,  we  show 
there  is  a  significant  increase  in  the  branching  of  the  self-assembled  network.  We 
also  show  that  magnetic  spheres  of  equal  magnetic  field  strength  can  segregate  due 
only  to  differences  in  their  field  shape.  We  find  that  the  segregation  increases  ap¬ 
proximately  linearly  with  acceleration  over  the  acceleration  range  studied.  We  also 
present  evidence  that  systems  of  vertically  vibrated  magnetic  spheres  of  equal  mass, 
size,  and  surface  properties  can  segregate  clue  only  to  differences  in  the  particles’ 
magnetic  field  strengths.  The  observed  segregation  is  shown  to  increase  with  the 
proportion  of  weaker  particles,  and  it  also  increases  approximately  linearly  with  ac¬ 
celeration  over  the  acceleration  range  studied.  Finally,  we  show  that  the  segregation 
is  accompanied  by  a  significant  decrease  in  magnetic  energy. 

Segregation  in  binary  mixtures  of  excited  granular  material  is  a  widely  ob¬ 
served  phenomenon  [59].  Segregation  due  to  differences  in  size  [39,  40,  60],  shape  [61, 


62],  density,  and  particle  surface  properties  [39]  has  been  observed  to  occur  under 
various  circumstances.  In  particular,  there  is  significant  current  interest  in  the  seg¬ 
regation  and  mixing  of  wet  granular  mixtures  as  an  example  of  segregation  among 
mutually  attractive  particles  [39,  40].  In  granular  systems,  the  addition  of  attractive 
forces  between  particles  has  been  observed  to  both  cause  and  prevent  segregation  un¬ 
der  different  circumstances  [39] .  Recent  observations  have  also  been  made  on  pattern 
formation  in  systems  of  dry  magnetic  spheres,  as  an  investigation  of  an  anisotropi- 
cally  attractive  granular  material  [26,  41].  Blair  and  Kudrolli  [26]  found  that,  in  a 
vertically  vibrated  2D  mixture  of  magnetic  and  non-magnetic  granular  spheres,  the 
magnetic  spheres  can  self  organize  and  cluster,  depending  on  the  volume  fraction 
of  magnetic  particles  and  on  the  relative  strength  of  the  magnetic  dipole-dipole  in¬ 
teraction  to  the  external  vibration  amplitude.  In  our  earlier  work  [41],  we  showed 
that  the  shape  of  the  magnetic  field  also  influences  into  what  pattern  the  beads 
sclf-assemble. 

In  addition  to  the  previous  work  on  vertically  vibrated  magnetic  granular 
spheres  [26,  41],  several  groups  have  investigated  vertically  vibrated  monolayers  of 
unmagnetized  granular  particles  [63,  64,  65].  However,  the  behavior  of  binary  sys¬ 
tems  of  granular  matter  that  differs  only  in  the  shape  or  strength  of  the  interaction 
potential,  to  our  knowledge,  has  not  been  investigated.  The  results  described  in 
this  chapter  can  therefore  either  be  viewed  as  an  investigation  of  a  novel  binary 
granular  material  with  anisotropic  interactions,  or  as  a  simple  model  system  of  a 
binary  dipolar  fluid. 
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5.2  Experimental  Details 


Figure  5.1  is  a  schematic  of  the  magnetic  particles  used  for  this  chapter.  We  investi¬ 
gate  particles  encasing  magnetic  cores  of  three  types  which  we  call  ‘weak’,  ‘strong’, 
and  ‘long’  magnets.  The  properties  of  all  particles  are  shown  in  table  5.1.  The  weak 
particles  have  nearly  identical  shape,  size,  and  mass  as  the  strong  particles,  with  the 
only  relevant  differing  property  being  their  relative  magnetic  strengths  (measured 
as  the  maximum  magnetic  field  on  the  surface  of  an  isolated  sphere).  The  strong 
and  lo7ig  particles  have  comparable  magnetic  strengths  and  are  nearly  identical  in 
shape,  size,  and  mass,  with  the  only  relevant  differing  property  being  their  mag¬ 
netic  core  shapes,  i.e.  their  magnetic  field,  shapes.  The  strong  and  weak  magnets 
were  encapsulated  with  non-magnetic  filler  of  appropriate  size  and  mass  so  that  the 
magnets  were  centrally  located  and  the  finished  particles  were  identical  in  mass. 
For  each  strong  and  weak  particle,  the  non-magnetic  filler  consisted  of  two  pieces 
of  hollow  cylindrical  brass  tubing  with  inner  and  outer  diameters  and  lengths  such 
that  the  magnets  were  rigidly  fixed  and  centrally  located,  and  the  final  masses  were 
well  matched  to  those  of  the  long  particles. 

The  container  was  constructed  with  a  black  Dclrin  bottom  and  a  mixture  of 
two  types  of  particles  was  placed  in  the  container  so  that  the  particles  were  well 
mixed,  i.e.  there  were  relatively  small  numbers  of  like-particle  contacts  prior  to 
vibration,  a  was  varied  from  4 g  to  10 g  and  the  system  was  imaged  from  above  at  a 
rate  of  1  Hz  with  a  1280x1024  pixel  color  camera. 
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Tab.  5.1:  A  table  of  particle  properties. 


Particle  Name 

Colors  (North/South) 

Mass  (g) 

Strength  (kG)“ 

Pb 

Z 

O  ZZZ 

l  (cm) 

Magnet  Material 

Weak 

White/Yellow 

5.4  ±0.1 

0.22  ±0.03 

1.0  ±0.1 

1.0  ±0.1 

0.52  ±0.1 

Ceramic 

Strong 

Bine /Green 

5.4  ±0.1 

0.60  ±0.08 

2.7  ±0.4 

2.6  ±0.4 

0.51  ±0.1 

NdFeB 

Long 

Yellow/ Green  or  Yellow 

5.3  ±0.1 

0.54  ±0.09 

2.2  ±0.3 

17  ±3 

1.42  ±0.1 

Ceramic 

a  Measured  as  the  maximum  surface  field  on  an  isolated  particle. 


6  Where  the  z  axis  is  parallel  to  the  line  connecting  the  poles,  with  the  origin  at  the  particle  center.  Moments  are  calculated  using  the  magnetic 


strength  and  the  core  shape,  and  are  normalized  by  moments  of  a  Weak  particle. 
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Fig.  5.1:  A  schematic  of  the  magnetic  particles  used,  together  with  a  qualitative  illustra¬ 
tion  of  their  respective  magnetic  field  lines  (obtained  using  Vizimag  software). 


5.3  Results 

Figure  5.2  is  a  representative  image  a  system  of  2  strong  and  55  long  particles  during 
vibration.  Note  that  3  to  4-fold  branching  occurs  around  the  short  particles,  while 
branching  is  uncommon  amongst  the  long  particles.  The  presence  of  a  relatively 
small  number  of  strong  particles  significantly  increases  the  branching  in  the  network 
pattern. 

In  systems  of  larger  mixture  fractions,  a  different  pattern  is  manifest.  As  an 
example  of  the  observed  patterns,  Fig.  5.3(a)  shows  the  mixed  initial  state  of  88 
weak  and  23  strong  particles.  The  initial  state  was  made  to  have  a  low  number  of 
like-particle  contacts.  Figure  5.3(b)  is  the  final  state  of  the  system  in  Fig.  5.3(a) 
after  being  shaken  at  9.0 g  for  300  seconds.  The  final  configuration  (Fig.  5.3(b)) 
is  highly  segregated,  with  all  strong  blue/green  (dark/dark  in  grayscale)  particles 
in  two  dense  clumps,  and  the  weak  white/yellow  (light/light  in  grayscale)  particles 
forming  a  loose,  branching  network  pattern.  The  evolution  of  this  system  as  shown 
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Fig.  5.2:  A  representative  image  of  a  system  of  2  strong  (blue/green)  and  55  long  (solid 
yellow)  particles  during  vibration  at  8.6 g.  Note  that  3  to  4-fold  branching  oc¬ 
curs  around  the  short  particles,  while  branching  is  uncommon  amongst  the  long 
particles. 
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in  Figs.  5.3(a)  and  (b)  indicates  that  systems  of  otherwise  identical  magnetic  spheres 
can  segregate  due  solely  to  differences  in  magnetic  held  strength. 

A  mixed  initial  state  of  51  strong  blue/green,  and  59  long  yellow/green  (light/dark 
in  grayscale)  particles  is  shown  in  Fig.  5.4(a).  After  300  seconds  of  excitation  at 
9.6 g,  the  system  evolves  into  the  configuration  shown  in  Fig.  5.4(b).  The  final  state 
(Fig.  5.4(b))  is  well  segregated,  with  the  strong  particles  having  a  tendency  toward 
hexagonal-close-packed  (HCP)  clumps,  and  the  long  particles  in  long  unbranched 
chains.  The  strong  particles  roughly  tend  toward  HCP  clumps  while  the  long  parti¬ 
cles  form  long,  unbranching  chains  as  a  consequence  of  their  differences  in  magnetic 
held  shape,  as  has  been  remarked  upon  earlier  in  a  study  of  pattern  formation  in 
monodisperse  2D  systems  of  magnetic  spheres  [41].  As  evidenced  in  Figs.  5.4(a) 
and  (b),  systems  of  otherwise  identical  magnetic  spheres  can  segregate  due  solely  to 
differences  in  magnetic  held  shape. 

For  further  analysis,  we  track  the  positions,  3D  dipole  orientations,  and  particle 
types  of  all  particles  during  our  experiments.  To  quantify  the  segregation  in  each 
system,  we  designed  a  program  that  finds  all  particle-particle  contacts  in  every 
image.  The  program  identifies  particles  as  ‘in  contact’  if  the  centers  of  the  particles 
are  within  9/8  particle  diameters  of  each  other.  We  tested  different  cutoff  values 
and  found  that  our  results  do  not  qualitatively  depend  on  the  cutoff.  We  then  define 
a  segregation  parameter  for  each  image, 

o  _  Caa  Cbb 

~  CAA  +  CAB  +  CBB  +  CAJ 

where  Ct]  is  the  total  number  of  ‘contacts’  between  particles  of  type  i  and  particles  of 


(5.1) 
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Fig.  5.3:  (a)  The  mixed  initial  state  of  a  system  of  88  weak  (white/yellow,  light/light  in 
grayscale)  and  23  strong  (blue/green,  dark/dark  in  grayscale)  particles,  (b)  The 
final  segregated  state  of  the  88  weak  and  23  strong  particles  after  being  shaken 
at  9.0 g  for  300  seconds.  Note  that  there  are  two  dense  clumps  of  the  strong 
particles,  while  the  weak  particles  form  a  loose  network  pattern. 
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Fig.  5.4:  (a)  The  mixed  initial  state  of  a  system  of  51  strong  (blue/green,  dark/dark  in 
grayscale)  and  59  long  (yellow/green,  light/dark  in  grayscale)  particles,  (b)  The 
final  state  of  the  51  strong  and  59  weak  particles  after  being  shaken  at  9.6 g  for 
300  seconds.  Note  that  the  strong  particles  are  largely  in  clusters  and  segregated 
from  the  long  particles,  which  form  long  unbranching  chains. 
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type  j  in  the  image  ( A  and  B  are  particle  type  labels).  This  segregation  measure  has 
the  advantages  of  symmetry  under  interchange  of  particle  type  and  equal  weighting 
of  the  segregation  of  each  particle  type.  0  <  S  <  2,  and  S  —  1,  on  average,  for 
random  mixtures  of  all  mixture  ratios. 

That  S'  =  1  for  random  mixtures  of  all  mixture  ratios,  can  be  seen  by  exam¬ 
ining  the  probabilities  of  contacts  of  each  type: 


Ni 


Pij  = 


Ni+N. 


Ni- 1 


Ni+Nj- 1’  1  3 

where  pij  is  the  probability  that  a  contact  associated  with  a  particle  of  type  i  is  with 
a  particle  of  type  j,  and  is  the  total  number  of  particles  of  type  i  in  a  system. 


So,  on  average, 


Cij  nNiPij , 


(5.2) 


where  k  is  the  mean  number  of  contacts  per  particle;  hence,  S  =  1. 

Figure  5.5(a)  is  a  plot  of  S  versus  time  for  a  single  trial  of  the  88  weak/23 
strong  system  shown  in  Fig.  5.3.  S  increases  from  an  initial  value,  0.86,  to  a  final 
value,  1.48,  indicating  that  the  system  began  in  a  state  which  was  ‘ supermixed,’’ 
(since  we  intentionally  placed  unlike  beads  next  to  each  other),  i.e.,  more  mixed 
than  a  random  mixture  ( S  <  1),  and  evolved  into  a  segregated  state  ( S  >  1). 
Figure  5.5(b)  is  a  plot  of  S  versus  time  for  a  single  trial  of  the  51  strong/59  long 
system  in  Fig.  5.4.  S  increases  from  an  initial  value,  0.88,  to  a  final  value,  1.39,  again 
indicating  that  the  system  evolved  from  a  mixed  state  to  a  segregated  state.  While 
trials  with  similar  particle  mixtures  and  acceleration  amplitudes  but  differing  initial 
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conditions  showed  similar  qualitative  long  time  segregation  behavior,  the  short  time 
behavior  of  S  versus  time  varied  significantly  in  different  trials. 

To  compare  the  degree  of  segregation  of  different  systems  after  300  seconds,  we 
define  a  parameter,  Sf,  as  the  mean  of  S  over  the  final  20  seconds  of  a  300  second- 
long  experiment.  Figure  5.6  (a)  is  a  plot  of  Sf  versus  the  strong  particle  fraction 
for  strong/weak  systems  of  111  particles  driven  at  9.0 g  (error  bars  are  shown  as  one 
standard  deviation  of  S,  calculated  over  the  final  20  seconds  of  a  single  experiment). 
Sf  decreases  as  the  fraction  of  strong  particles  is  increased.  Figure  5.6  (b)  is  a  plot 
of  Sf  versus  acceleration,  a/g ,  for  a  system  of  67  weak  and  44  strong  particles  (filled 
circles)  and  a  system  of  67  strong  and  44  long  particles  (filled  triangles).  Sf  increases 
roughly  linearly  with  acceleration  for  both  systems.  Note  also,  there  are  several  data 
points  in  Fig.  5.6  (b),  where  Sf  <  1,  indicating  that,  after  300  seconds,  the  systems 
were  supermixed.  This  is  not  an  indication  that  the  systems  prefer  a  supermixed 
state;  rather,  it  is  due  both  to  the  fact  that  the  systems  began  in  a  supermixed 
state  and  to  the  particles’  lack  of  mobility.  The  particles’  lack  of  mobility  at  low 
accelerations  is  likely  because  the  average  input  (driving)  energy  is  less  than  the 
magnetic  particle-particle  interaction  energy.  In  no  case,  for  any  of  our  systems, 
have  we  observed  a  significant,  non-transient  decrease  in  S  in  time. 

5.4  Energetic  Considerations 

To  gain  insight  into  the  cause  of  the  observed  segregation,  we  numerically  examined 
the  magnetic  particle  model  described  in  Chapter  2.  We  vary  q  with  h  so  that  the 
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Fig.  5.5:  a)  The  segregation,  S,  versus  time  for  the  system  of  88  weak  and  23  strong 
particles  shown  in  Fig.  5.3.  Note  that  S  increases  significantly  over  time,  and  it  is 
greater  than  1  (the  S'  of  a  randomly  mixed  system)  at  the  end  of  the  experiment, 
b)  S  versus  time  for  the  system  of  51  strong  and  59  long  particles  shown  in 
Fig.  5.4.  Note  that  S  again  increases  significantly  over  time,  and  it  is  greater 
than  1  at  the  end  of  the  experiment 
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Fig.  5.6:  (a)  Sf  versus  the  strong  particle  mixture  fraction  for  strong/weak  systems  of  111 
particles  accelerated  at  9.0 g.  Note  that  Sf  decreases  as  the  number  of  strong 
particles  in  the  mixture  increases,  (b)  Sf  versus  acceleration,  a/g,  for  systems 
of  111  particles.  Mixtures  of  44  strong  and  67  weak  particles  are  shown  as 
filled  circles,  while  mixtures  of  67  strong  and  44  long  particles  are  shown  as 
filled  triangles.  Note  that,  for  both  mixtures,  Sf  increases  roughly  linearly  with 
acceleration. 
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dipole  moment,  2 qh,  is  constant.  This  model  is  a  generalized  version  of  the  two- 
charge  model  used  previously  [41,  27];  it  reduces  to  the  two-charge  model  in  the 
limit  of  w/D  — >  0.  It  can  also  be  reduced  to  a  two-point-dipole  model  in  the  limit 
of  of  h/D  — »  0,  and  a  single,  centrally  located  point  dipole  model  if  h/D  0  and 
w/D  — >  0. 

By  considering  our  cylindrical  magnets  as  two  uniformly  charged  discs  of  oppo¬ 
site  sign,  but  equal  total  charge  magnitude,  we  can  then  attempt  to  match  magnetic 
scalar  potentials  of  our  model  particles  to  our  real  particles.  Using  our  simple  two- 
parameter  model,  we  match  the  magnetic  scalar  potential  of  our  real  and  model 
particles  up  to  the  octopolc  terms.  The  weak  and  strong  particles,  having  similar 
shaped  cores,  are  both  modeled  with  two  parallel  point  dipoles  separated  by  a  dis¬ 
tance,  w  =  0.31.D,  as  shown  in  Fig.  5.7(b).  The  point  dipoles  are  approximated 
using  h  =  0.01.D  for  both  the  weak  and  strong  particles,  with  q  =  11.5  and  q  =  34.5, 
respectively,  to  account  for  the  difference  in  magnetic  strength.  The  long  parti¬ 
cles  are  modeled  with  two  opposite  point  charges  of  magnitude  2.0,  separated  by  a 
distance  of  0.69 D,  as  shown  in  Fig.  5.7(c). 

Using  our  experimental  particle  type,  position,  and  orientation  data  extracted 
from  the  images,  we  use  our  model  to  approximate  the  total  magnetostatic  energy 
of  our  systems  for  every  frame.  We  use  the  usual  Coulomb  energy  formula,  Uij  = 
,  and  sum  over  all  interparticle  pairs  of  charges  (excluding  the  energy  of  pairs 
of  charges  within  the  same  particle)  to  obtain  the  energy  of  our  experimentally 
observed  patterns  of  particles,  Uex.  As  a  reference,  we  use  the  magnetic  energy  of 
a  randomly  mixed  state  of  the  same  spatial  and  orientational  structure,  which  we 
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Fig.  5.7:  (a)  A  schematic  of  our  general  magnetic  particle  model,  as  introduced  in  Fig.  2.3. 

(b)  For  the  weak  and  strong  particles,  the  model  reduces  to  a  two  point-dipole 
model,  with  the  dipoles  separated  by  a  distance,  w  =  0.31.D,  where  D  is  the 
particle  diameter,  (c)  For  the  long  particles,  the  model  in  reduces  to  a  two- 
charge  model,  with  the  charges  separated  by  a  distance,  h  =  0.69 D. 

obtain  by  taking  the  particle  tracking  data  from  the  experiment,  and  interchanging 
the  particle  types  of  randomly  chosen  particle  pairs  (keeping  the  particle  locations 
and  dipole  orientations  fixed).  The  magnetic  energy  was  then  recalculated.  This 
randomization  and  calculation  process  was  repeated  eleven  times  and  the  mean 
magnetic  energy  of  the  randomized  system,  Ur  was  calculated  for  each  frame.  After 
eleven  randomizations,  the  time  average  of  the  standard  deviation  of  the  randomized 
magnetic  energy  was  5%  of  Ur. 

In  Fig.  5.8  (a),  we  have  plotted  Uex  (filled  circles)  versus  time  for  the  system  of 
88  weak  and  23  strong  particles  shown  in  Fig.  5.3.  Ur  (open  circles)  is  also  plotted 
in  Fig.  5.8  (a).  The  standard  deviation  of  Ur  is  essentially  constant  in  time  at  ±5, 
approximately.  Note  that  Uex  decreases  over  time,  while  Ur  is  roughly  constant;  also 
Uex  <  Ur  for  all  time.  The  noticeable  decrease  of  Uex  below  Ur  over  time  indicates 
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that  the  decrease  in  Uex  is  not  an  artifact  of  indiscriminant  pattern  formation,  but 
rather  a  part  of  the  segregation  phenomenon.  Fig.  5.8  (b)  is  a  plot  of  Uex  —  Ur 
versus  time.  Note  that  Uex  —  Ur  decreases  steadily  after  the  first  few  seconds. 
The  outlying  points  in  Fig.  5.8  (a)  and  (b)  are  due  to  errors  in  particle  position 
and  orientation  extraction,  where  false  particles  were  detected  or  orientations  were 
improperly  extracted.  While  the  observed  segregation  patterns  tend  toward  states 
of  low  magnetostatic  energy,  it  is  not  expected  a  priori  that  a  driven  dissipative 
system  would  arrange  into  a  globally  minimal  magnetostatic  energy  state  over  long 
times.  Rather,  the  observed  patterns  can  be  understood  as  non-equilibrium  steady 
states. 


5.5  Conclusions 

Using  a  simple  experiment  of  a  model  system,  we  have  shown  that  in  dilute  self- 
assembled  networks,  a  small  change  in  the  multipole  moment  distribution  of  a  small 
fraction  of  the  particles  can  significantly  increase  the  branching  in  the  self-assembled 
network.  This  has  significant  implications  for  the  macroscopic  properties  of  polymer 
networks,  and  it  gives  insight  into  the  possible  properties  of  branching  particles  in 
biologically  relevant  systems.  Our  results  also  show  that  otherwise  identical  sys¬ 
tems  of  magnetized  granular  particles  can  segregate  due  only  to  differences  in  their 
magnetic  fields.  Vertically  vibrated  monolayers  of  magnetic  spheres  can  segregate 
both  by  held  strength  and  by  held  shape.  In  binary  systems  of  particles  with  differ¬ 
ing  held  strength,  and  all  other  properties  identical,  the  segregation  increases  with 
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Fig.  5.8:  a)  The  magnetic  energy,  Uex,  versus  time  for  the  system  of  88  weak  and  23  strong 
particles  shown  in  Fig.  5.3  is  plotted  as  filled  circles.  The  magnetic  energy  of  the 
randomly  mixed  state,  Ur,  is  also  plotted  as  open  circles  for  reference.  Note  that 
Uex  decreases  quickly  within  the  first  ten  seconds,  then  decreases  more  gradually 
over  time,  while  Ur  remains  roughly  constant,  b)  A  plot  of  Uex  —  Ur  versus  time. 
Uex  —  Ur  decreases  steadily  after  the  first  few  seconds. 


* 


71 


the  proportion  of  weaker  particles,  and  it  also  increases  approximately  linearly  with 
acceleration  over  the  acceleration  range  studied.  Binary  systems  of  particles  with 
differing  field  shape  and  all  other  properties  identical  also  show  approximately  linear 
increases  in  segregation  as  acceleration  is  increased  over  the  acceleration  range  stud¬ 
ied.  Segregation  occurs  in  conjunction  with  a  decrease  in  magnetic  energy,  with  the 
energy  decrease  being  mostly  due  to  the  actual  segregation,  rather  than  an  evolution 
of  the  spatial  pattern. 


5.6  Outlook 

The  weak  plastic  construction  of  our  particles  prevents  us  from  probing  the  behavior 
of  our  systems  at  higher  accelerations;  however,  it  seems  likely  that  at  high  accel¬ 
erations,  the  input  vibration  energy  would  exceed  the  magnetic  interaction  energy, 
and  systems  such  as  ours  would  remix.  While  we  expect  that  the  qualitative  behav¬ 
ior  of  larger  systems  to  be  similar,  the  small  size  of  the  systems  we  studied  makes 
it  difficult  to  obtain  high-quality  quantitative  statistical  results.  Possible  future 
work  includes  studies  of  larger  binary  systems  of  magnetic  particles,  and  using  more 
sturdily  constructed  particles  that  allow  the  investigation  of  higher  accelerations. 

It  would  also  be  interesting  to  investigate  the  interaction  properties  of  parti¬ 
cles  which  induce  branching  in  biological  self-assembly  and  polymer  networks,  and 
relate  their  properties  to  those  of  our  magnetic  branching  particles.  The  interac¬ 
tions  between  biological  macromolecules,  and  methods  of  quantitative  comparison 
of  interactions  will  be  discussed  in  Chapter  6.  It  may  be  possible  to  simulate  known 
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branching  particles  and  other  biological  macromolecnles  using  methods  described 
in  Chapter  6,  and  it  would  also  be  intersting  to  insert  new,  synthetic  particles  in 
the  simulation  to  investigate  their  roles  in  the  self-assembly  of  biologically  relevent 
proteins.  It  is  possible  that  insight  gained  from  such  exploration  could  give  insight 
into  possible  real  synthetic  particles  that  would  have  similar  effects  on  real  biological 
self-assembly. 
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6.  PROTEINS:  ANISOTROPICALLY  INTERACTING 


SELF-ASSEMBLING  PARTICLES 

6.1  Self-Assembly  in  Biological  Systems 

Self-assembly  is  of  particular  interest  in  molecular  biology  [36,  37].  A  wide  array  of 
biologically  relevant  systems  form  through  self-assembly  [66].  Protein  self-assembly 
is  of  particular  interest  because  of  the  many  diseases  that  are  a  result  of  the  self- 
assembly  [1,  2,  3,  4,  5,  36,  67].  Hemoglobin  S,  tubulin,  and  collagen  self-assemble  in 
sickle  cell  anemia,  cancerous  tumor  growth,  and  brittle  bone  disease,  respectively. 
Hemoglobin  S,  a  mutation  of  the  normal  hemoglobin  A  protein,  self- assembles  into 
long  twisting  rope-like  structures  in  the  red  blood  cells  of  diseased  individuals,  and 
inhibits  the  transportation  of  oxygen  within  the  body  [36]. 

6.1.1  Hemoglobin 

Hemoglobin  is  a  protein  mainly  found  in  red  blood  cells;  its  chief  function  is  the 
transport  of  oxygen  within  the  body.  Hemoglobin  A  is  the  normal  form  of  the 
protein,  while  hemoglobin  S  is  a  variant  which  results  from  a  small  mutation  of  the 
protein  structure.  While  hemoglobin  A  typically  does  not  self-organize,  hemoglobin 
S  forms  long  twisting  rope-like  structures  inside  red  blood  cells  [36].  Figure  6.1(a) 


is  a  rendering  of  the  hemoglobin  A  protein  taken  from  the  known  crystallographic 
structure  [68],  while  fig.  6.1(b)  is  a  rendering  of  the  hemoglobin  S  protein  [69].  Note 
that,  in  fig.  6.1,  hemoglobin  S  appears  to  be  two  hemoglobin  A  proteins  bound 
together;  the  binding  of  the  two  hemoglobin  A-like  proteins  is  a  consequence  of 
the  mutation.  Figure  6.2  is  a  schematic  of  a  typical  self-assembly  of  hemoglobin  S 
proteins.  The  self-assembly  of  these  proteins  is  a  direct  result  of  the  alteration  of 
the  interactions  from  the  normal  to  mutated  case. 

6.1.2  Tubulin 

Tubulin  is  a  roughly  cylindrical  dimeric  protein  which  self-assembles  into  long,  hol¬ 
low,  cylindrical  structures  called  microtubules  inside  most  living  cells  [70,  71].  Fig¬ 
ure  6.3(a)  is  a  rendering  of  the  tubulin  protein,  drawn  from  the  known  structure  of 
the  protein  from  crystallography  [72];  fig.  6.3(b)  is  a  schematic  of  a  normal  human 
microtubule.  Microtubules  play  a  crucial  role  in  the  transport  of  materials  within 
the  cell,  cell  structure,  and  cell  division.  Because  of  the  key  role  microtubule  growth 
has  in  cell  division,  controlling  tubulin  self-assembly  is  the  purpose  of  many  modern 
chemotherapy  drugs  [2], 


6.1.3  Collagen 

Collagen  is  a  protein  which  plays  and  essential  role  in  the  formation  of  bones  and 
cartilage  in  the  body.  There  are  many  variants  of  the  protein,  some  of  which  are 
mutants  which  cause  bone  and  joint  disease  [67].  Figure  6.4(a)  is  a  rendering  of  a  col¬ 
lagen  protein,  drawn  from  its  known  structure  from  crystallography  [73];  fig.  6.4(b) 
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Fig.  6.1:  a)  A  rendering  of  the  electrostatic  surface  potential  of  the  hemoglobin  A  pro¬ 
tein  drawn  from  the  known  crystallographic  structure  using  Protein  Explorer 
Software.  Red  indicates  an  area  of  relative  negative  charge,  blue  is  an  area  of 
positive  charge,  and  white  is  neutral,  b)  A  rendering  of  the  hemoglobin  S  protein 
by  similar  method.  Note  that  hemoglobin  S  appears  to  be  two  hemoglobin  A 
proteins  bound  together. 
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Frank  Ferrone,  Drexel  University 


Fig.  6.2:  A  schematic  of  the  twisted  rope-like  structure  of  self-assembled  hemoglobin  S. 
Color  has  no  intended  significance.  Schematic  drawn  by  Frank  Ferrone. 

is  a  confocal  fluorescence  microscopy  image  of  the  self-assembled  network  structure 
formed  by  normal  collagen  proteins  [74],  It  is  believed  that  the  mutated  collagen 
associated  with  disease  causes  a  significant  change  in  collagen  interactions,  which 
alters  the  microscopic  and  macroscopic  physical  properties  of  the  self-assembled 
network  [67].  Understanding  the  roles  of  the  protein  interactions  on  the  pattern  of 
self-assembly  is  of  significant  interest  because  of  the  possible  value  for  therapeutics. 

6.2  Protein  Electrostatics 

Protein  interactions,  while  anisotropic  in  nature,  are  generally  far  too  complicated 
to  be  described  by  a  point  dipole.  Interactions  between  proteins  can  be  divided, 
principally,  into  two  types  [75]: 

1.  Hydrophobic  interactions  due  to  the  relative  affinities  to  water  of  the  different 
surfaces  of  the  protein  [76]. 


77 


Fig.  6.3:  a)  A  rendering  of  the  electrostatic  surface  potential  of  tubulin,  drawn  from  the 
known  crystallographic  structure  using  Protein  Explorer  Software.  Red  indicates 
an  area  of  relative  negative  charge,  blue  is  an  area  of  positive  charge,  and  white  is 
neutral,  b)  A  schematic  of  a  normal  human  self-assembled  structure  of  tubulin, 
a  microtubule. 
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Fig.  6.4:  a)  A  rendering  of  the  electrostatic  surface  potential  of  collagen,  drawn  from  the 
known  crystallographic  structure  using  Protein  Explorer  Software.  Red  indicates 
an  area  of  relative  negative  charge,  blue  is  an  area  of  positive  charge,  and  white 
is  neutral,  b)  A  confocal  fluorescence  microscopy  image  of  the  self-assembled 
network  structure  formed  by  normal  collagen.  Image  taken  by  Sejin  Han  and 
Wolfgang  Losert 
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2.  Electrostatic  interactions  due  to  partial  electron  charges  that  are  induced  in 
the  protein  by  the  surrounding  solution. 

Although  both  types  of  interactions  are  important  for  most  self-assembly,  the  elec¬ 
trostatic  properties  of  the  proteins  are  more  easily  quantified,  and  thus,  I  will  focus 
upon  the  protein  electrostatics  for  the  remainder  of  this  chapter. 

Using  X-ray,  electron  crystallography,  and  nuclear  magnetic  resonance,  atom¬ 
ically  detailed  structure  information  for  many  proteins  has  been  detrermined  and  is 
availible  online  at  the  Protein  Data  Bank  (PDB)  website,  http://www.rcsb.org/pdb/. 
The  PDB  contains  separate  hies  (.pdb  hies)  for  each  protein,  with  each  hie  providing 
the  positions  of  the  atoms  in  a  protein.  Another  online  resource,  the  PDB2PQR 
Server  at  http://agave.wustl.edu/pdb2pqr/,  takes  the  .pdb  hies  as  input  and  uses 
molecular  modeling  (AMBER,  CHARMM,  or  PARSE)  to  determine  approximate 
partial  charges  of  each  atom  in  the  protein,  assuming  the  protein  is  under  normal 
physiological  conditions.  The  PDB2PQR  Server  outputs  a  hie  (a  .pqr  hie)  which 
contains  the  position  and  charge  of  each  atom  in  the  .pdb  hie.  Using  the  atomic 
position  and  charge  data,  electrostatic  moments  of  many  proteins  can  be  calculated. 
I  have  calculated  the  electrostatic  moments  of  hemoglobin  A  (protein  database,  pdb, 
code  1A3N),  hemoglobin  S  (pdb  code  2HBS),  tubulin  (pdb  code  1JFF),  and  collagen 
(pdb  code  1CAG);  the  independent  moments  of  these  proteins,  up  to  quadrupole 
terms,  are  shown  in  Table  6.1,  and  the  octopole  moments  are  shown  in  Table  6.2. 
Note  that  the  dipolar  and  higher  order  moments  of  all  proteins  listed  in  Tables  6.1 
and  6.2  are  all  non-zero,  indicating  that  the  proteins  interact  anisotropically.  In 
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particular,  the  total  dipole  moment  of  tubulin  is  293 e~nm,  or  1410  Debye;  this 
implies  that,  neglecting  the  effects  of  the  surrounding  fluid,  the  dipole-dipole  inter¬ 
action  energy  between  tubulin  proteins  is  on  the  order  of  100 kT  when  the  proteins 
are  in  close  proximity.  Our  calculated  charge  and  dipole  moments  of  tubulin  are 
similar  to  those  calculated  by  J.  A.  Brown  [77],  while  our  calculated  dipole  moments 
of  hemoglobin  A  is  significantly  lower  than  the  480D  dipole  moment  calculated  by 
S.  Takashima  [78].  The  discrepancy  between  our  hemoglobin  A  dipole  moment  and 
that  calculated  by  Takashima  is  likely  due  to  our  use  of  the  higher  resolution  crystal¬ 
lographic  structure  of  the  protein  and  also  due  to  our  use  of  the  improved  molecular 
modeling  package  (AMBER)  implicit  in  the  PDB2PQR  Server.  The  use  of  other 
available  molecular  modeling  packages  (CHARMM  and  PARSE)  in  the  PDB2PQR 
Server  did  not  significantly  change  the  calculated  moments. 

6.3  Multipolar  Models 

Although  the  electrostatic  potentials  of  proteins  are  somewhat  simple  to  estimate 
using  the  atomically  detailed  position  and  charge  data,  because  of  the  large  num¬ 
ber  of  partial  charges  associated  with  the  proteins  (on  the  order  of  104,  typically), 
the  computation  of  electrostatic  interactions  between  proteins  can  be  computation¬ 
ally  expensive.  The  relation  between  the  computation  time,  t,  of  forces  between 
electrostatically  interacting  particles  and  the  number  of  charges,  N,  is  t  oc  N2, 
if  no  cutoff  length  is  introduced.  Because  of  the  large  number  of  charges  associ¬ 
ated  with  most  proteins,  and  the  computational  time  required  to  properly  calculate 
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Tab.  6.1:  Charge,  dipole,  and  quadrupole  electrostatic  moments  of  hemoglobin  A,  hemoglobin  S,  tubulin,  and  collagen.  Ptotai  is  in  units 
of  Debyes  and  all  other  moments  are  calculated  in  units  of  electron  charge  and  nanometers. 


Protein 

pdb  code 

Q 

Ptotai  (Debyes) 

Px 

Py 

Pz 

Qxx 

Qyy 

Qxy 

Qxz 

Qyz 

Hemoglobin  A 

1A3N 

11.8 

215 

1.12 

-3.22 

-2.91 

13.1 

-58.7 

-2.39 

-27.7 

5.16 

Hemoglobin  S 

2HBS 

28.8 

545 

8.91 

3.92 

-5.83 

-11.7 

209 

327 

85.4 

100 

Tubulin 

1JFF 

-18 

1410 

7.43 

-27.5 

7.07 

-26.6 

48.4 

154 

105 

84.9 

Collagen 

1CAG 

0 

1150 

-24.0 

-.06 

-.04 

-15.9 

8.38 

-1.42 

1.48 

-.072 

Tab.  6.2:  Octopole  electrostatic  moments  of  hemoglobin  A,  hemoglobin  S,  tubulin,  and  collagen.  All  moments  are  calculated  in  units  of 
electron  charge  and  nanometers. 


Protein 

pdb  code 

0  XXX 

0Xxy 

0 XX  z 

Oxyy 

0Xyz 

O XZZ 

Oyyy 

Oyyz 

Oyzz 

O  zzz 

Hemoglobin  A 

1A3N 

-10.8 

11.7 

-14.9 

-10.8 

-17.4 

8.51 

-8.05 

-5.6 

-0.822 

-17.5 

Hemoglobin  S 

2HBS 

76.6 

92.5 

-16.0 

160 

-2.06 

-35.9 

110 

-22.0 

-17.8 

-42.6 

Tubulin 

1JFF 

83.3 

-140 

-44.0 

27.9 

78.0 

-103 

-146 

-0.753 

-47.0 

47.7 

Collagen 

1CAG 

-386 

-4.22 

-0.997 

-1.98 

-0.171 

-0.963 

0.348 

-0.0608 

-0.045 

-0.0721 

electrostatic  interactions,  large-scale  simulations  of  protein  self-assembly  based  on 
physical  protein-protein  interactions,  to  my  knowledge,  have  not  been  done.  Using 
electrostatic  moments  as  calculated  in  the  previous  section,  and  a  simple  Monte 
Carlo  algorithm,  I  am  able  to  develop  simple  8-charge  models  of  the  proteins  which 
match  the  electrostatic  moments  of  the  real  proteins  up  to  the  octopole.  This  is  a 
significant  step  toward  a  computationally  inexpensive  physical  interaction  model  of 
self- assembling  proteins. 

In  the  Monte  Carlo  algorithm,  8  charges  are  allowed  to  move  independently 
in  four-dimensional  position-charge  space  (q,  x,  y,  z ),  with  the  range  of  each  position 
parameter  (x,  y,  z)  constrained  such  that  the  charges  are  confined  within  the  model 
protein’s  spacial  boundaries,  q  is  also  constrained  such  that  \q\  <  qmax,  where 
Qmax  ~  100e_,  typically.  A  rough  outline  of  the  Monte  Carlo  algorithm  follows. 
Inputs: 

•  All  electrostatic  moments  of  the  protein  up  to  octopole. 

•  Range  of  (q,  x,  y,  z)  parameters. 

•  A  ‘thermal’  parameter,  kT.  Typically,  0.0001  <  kT  <  1. 

Internal  Parameters: 


The  variance, 


y=  E 

all  moments 


Enodel  -^real 


M 


real 


(6.1) 


where  Mrn0(.j(,]  and  Mrca]  are  electrostatic  moments  of  the  model  and  real 
protein,  respectively.  There  are  three  variance  parameters  in  the  algorithm, 
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^ new  i 


Voldi  Slid  Vmin • 


Vpld  —  Vnew 

•  The  Monte  Carlo  parameter,  T  =  e  kTVmin 


Steps: 

1.  V0id  and  Vmin  are  initialized  to  some  large  number  (typically  24). 

2.  All  q,  x,  y,  z  are  initialized  to  zero  (thus,  all  model  moments  are  also  initialized 
to  zero). 

3.  One  of  the  32  position-charge  parameters,  p,  is  chosen  at  random. 

4.  A  new  value,  pnew,  of  the  chosen  parameter  is  selected  from  a  uniform  distri¬ 
bution  with  the  range  as  input. 

5.  The  moments  of  the  model  are  calculated  and  Vnew  is  calculated. 

6.  A  random  number,  0,  between  0  and  1  is  chosen  from  a  uniform  distribution. 

7.  If  T  >  0,  Pnew  I®  kept  [p  Pnew)  Slid  Vo/J  Vnew  Also,  if  Vnew  Vnin  ■ 
Vmin  Vnew 

8.  The  program  loops  back  to  step  4. 

The  program  is  allowed  to  continue  to  loop  until  Vmtn  is  sufficiently  small.  Typ¬ 
ically,  runtime  is  about  24  hours  for  4  x  109  loops,  where  Vmin  ~  10~6  at  completion, 
such  that  all  moments  are  matched  to  better  than  1%.  Because  this  algorithm  takes 
only  electrostatic  moments,  the  ranges  of  parameters,  and  a  thermal  parameter  as 
input,  it  can  be  used  to  develop  minimal  electrostatic  models  of  nearly  any  protein. 
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The  algorithm  can  be  easily  be  modified  to  give  the  desired  model  accuracy.  The 
number  of  moments  to  match  determines  the  number  of  charges  needed  in  the  model 
(1  charge  for  monopole,  2  for  dipole,  4  for  quadrupole,  etc.). 

6.4  Outlook 

The  relatively  low  runtime  for  the  above  descried  algorithm  makes  it  especially 
useful  for  large-scale  self-assembly  simulations  where  electrostatic  calculations  are 
computationally  expensive.  However,  the  hydrophobic  interactions  must  also  be 
accounted  for  if  a  structure-based  model  of  protein  interactions  is  to  be  used.  If 
simple  models  of  the  hydrophobic  interactions  can  be  developed,  it  may  be  possible 
to  accurately  and  simply  model  total  protein  interactions  with  a  superposition  of 
minimal  electrostatic  and  hydrophobic  models.  One  possible  method  of  simulating 
hydrophobic  interactions  is  to  use  ‘charges’  that  have  Lennard-Jones  interactions. 
The  magnitudes  and  positions  of  such  ‘charges’  could  be  varied  and  possibly  matched 
with  information  extracted  from  the  Protein  Data  Bank,  similar  to  the  method 
used  in  the  algorithm  described  above.  However,  for  the  remainder  of  this  thesis  I 
investigate  fundamental  anisotropic  interaction  models  with  no  particular  reference 
to  hydrophobic  or  electrostatic  protein  modeling. 
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7.  THE  EFFECTS  OF  MULTIPOLAR  INTERACTIONS  ON 


SELF-ASSEMBLY 

7. 1  Overview 

As  illustrated  in  previous  chapters,  anisotropically  interacting  particles  form  a  wide 
variety  of  self-assembly  patterns.  Controlling  the  self-assembly  and  the  types  of 
patterns  that  are  formed  is  currently  a  major  goal  for  several  fields  of  study  ranging 
from  polymer  engineering  to  cellular  biology  and  medicine.  In  particular,  control¬ 
ling  protein  self-assembly  is  of  significant  interest  to  the  pharmaceutical  industry 
because  of  possible  therapeutic  value.  For  example,  the  chemotherapy  drug  Taxol 
introduces  additional  interactions  among  tubulin  proteins  in  cancer  cells  and  which, 
in  turn,  prevents  cell  division  and  the  growth  of  tumors  [2],  However,  because  of 
the  computational  expense  of  simulating  such  complex  self-assembly  systems  at  the 
molecular  level  and  a  desire  to  understand  the  self-assembly  at  a  fundamental  level, 
there  has  been  recent  renewed  interest  in  the  self-assembly  of  particles  with  simple 
interaction  potentials  such  as  dipolar  spheres  [29,  41]. 

In  our  recent  work  [41,  42],  we  investigated  the  effects  of  higher  order  multi¬ 
pole  moments  on  self-assembly  using  numerical  modeling  and  simple  experiments 
on  a  macroscopic  toy  model  of  magnetic  particles.  We  showed  that  a  small  change 


in  the  multipole  moment  distribution  of  the  particles  can  significantly  alter  the  pre¬ 
ferred  self-assembly  pattern  in  binary  and  monodisperse  systems.  In  this  chapter, 
we  show  that  if  the  interactions  of  all  particles  are  uniformly  altered,  complex  self- 
assembled  structures  can  be  stabilised.  As  an  example,  we  propose  a  simple  point 
charge  model  of  the  protein  tubulin,  and  we  show  that,  without  altering  the  dipole 
moment,  the  addition  of  higher  order  multipole  moments  to  the  model  can  ener¬ 
getically  stabilise  the  biologically  correct  self-assembly  pattern.  We  also  present  a 
physical  macroscopic  magnetic  analog  of  our  model,  and  show  that  it  has  many  of 
the  same  structural  characteristics  as  microtubules. 

7.2  Biological  Background 

As  illustrated  in  Chapter  6,  microtubules  are  long,  hollow,  cylindrical  structures 
formed  by  the  self-assembly  of  roughly  cylindrical  dimeric  tubulin  proteins  inside 
cells  [70,  71].  In  the  presence  of  different  drugs  and  ions,  tubulin  proteins  (each 
about  8.0  nm  long  and  4.6  nm  in  diameter)  self  assemble  into  a  plethora  of  other 
structures  including  rings,  sheets  [79],  spirals,  and  ‘macrotubules’  [37].  While  many 
variations  of  microtubule  structure  have  been  seen  [80],  for  in  vivo  human  micro¬ 
tubules,  tubulin  arranges  into  thirteen  long  filaments  ( protofilaments )  that  form  the 
wall  of  the  cylindrical  tube  structure.  In  the  microtubule  lattice,  protofilaments 
are  staggered  by  0.93  nm  from  neighboring  protofilaments  [81]  and,  at  the  edges  of 
depolymerizing  microtubules,  the  protofilaments  are  believed  to  fray  out  from  the 
main  tube  structure  and  break  off  [82], 


The  unusual  pattern  of  microtubule  growth  and  sudden  shrinkage  (termed  dy¬ 
namic  instability )  is  not  entirely  understood,  and  is  a  topic  of  active  research  [75,  83]. 
It  is  our  belief  that  simulations  of  tubulin  self-assembly  based  on  physical  interac¬ 
tions  between  proteins  could  elucidate  the  microtubule  growth  process,  and  perhaps 
aid  in  novel  drug  discovery.  However,  because  of  the  significant  computational  ex¬ 
pense  of  using  atomically  detailed  tubulin  interaction  models,  large  scale  simulations 
of  tubulin  self-assembly  based  on  physical  interactions,  to  our  knowledge,  have  not 
been  done.  The  goal  of  our  work  is  to  develop  an  interaction  model  which  accu¬ 
rately  reproduces  the  self-assembly  and  is  also  simple  enough  to  use  in  large-scale 
simulations. 

Many  investigations  have  been  done  in  an  attempt  to  understand  the  inter¬ 
action  between  individual  tubulin  dimers  and  its  relationship  to  microtubule  struc¬ 
ture  [75,  84,  85].  In  particular,  recent  computational  work  by  Sept  et  al.  [75]  correctly 
predicted  the  microtubule  filament  staggering,  and  gave  estimates  of  interfilament 
and  intrafilament  tubulin-tubulin  binding  energies  using  the  known  structure  of  mi¬ 
crotubules  [86].  However,  because  of  the  complexity  of  the  tubulin  protein  and 
the  computational  expense  of  accurately  reproducing  all  interactions,  the  study  was 
limited  to  static  energy  calculations  of  only  five  proteins  in  a  limited  number  of 
spacial  configurations. 

It  is  well  known  that  tubulin  is  a  protein  that  carries  a  strong  electrostatic 
dipole  moment  [77].  Recent  calculations,  based  on  the  known  structure  of  tubulin 
indicates  the  strength  of  the  intrinsic  dipole  moment  of  tubulin  is  about  1700D  [77]; 
neglecting  the  effects  of  the  surrounding  solution,  this  indicates  the  tubulin  dipole- 


89 


dipole  interaction  energy  is  on  the  order  of  100  kT.  It  has  also  been  shown  that  in 
vitro  microtubules  align  in  the  presence  of  an  electric  held  [87,  88] .  The  electrostatic 
polarity  of  microtubules,  and  its  possible  relationship  to  its  cellular  functionality  is 
a  topic  of  current  research  [89].  In  addition  to  the  known  effects  of  electric  fields  on 
microtubules,  it  has  also  been  shown  that  magnetic  [87]  and  even  gravitational  fields 
[90]  can  significantly  effect  microtubule  growth.  While  models  based  on  the  dipolar 
nature  of  tubulin  have  been  used  in  investigations  of  microtubule  function  [38,  77, 
91,  92],  models  which  investigate  the  role  of  higher  order  moments  in  tubulin  self- 
assembly  have  not  been  studied,  to  our  knowledge. 

7.3  Experimental  Toy  Model 

Figure  7.1  is  a  photograph  of  patterns  formed  by  cylindrical  permanent  magnets 
(dipolar  cylinders)  of  length  1.42  cm  and  diameter  0.95  cm,  similar  to  those  used  in 
the  ‘long’  magnetic  particles  in  Chapters  3-5.  The  magnets  form  a  staggered,  frayed 
yet  stable  microtubule- like  2D  structure  when  chains  of  the  magnets  are  brought 
together  with  consistent  orientation.  When  chains  of  magnets  are  brought  together 
with  opposite  orientation,  they  form  a  stable,  non- frayed  sheet,  much  like  the  zinc- 
stabilized  tubulin  sheets  seen  in  previous  work  [37].  Chains  of  magnets  can  also  be 
bent  to  form  a  closed  ring  shape,  much  like  the  observed  tubulin  rings  [37,  93].  In¬ 
spired  by  these  observations  of  self-assembled  dipolar  cylinders,  we  propose  a  simple 
interaction  model  that  reproduces  many  of  the  key  features  of  microtubules,  in¬ 
cluding  0.93  nm  staggering  between  neighboring  protofilaments,  stable  microtubule 
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Fig.  7.1:  A  photograph  of  patterns  formed  by  cylindrical  permanent  magnets  of  length 
1.42  cm  and  diameter  0.95  cm.  The  magnets  form  a  staggered,  frayed  yet  stable 
microtubule-like  structure  when  chains  of  the  magnets  are  brought  together  with 
consistent  orientation.  Chains  of  magnets  can  also  be  bent  to  form  a  closed  ring 
shape,  much  like  the  observed  tubulin  rings. 

structure,  and  fraying  protofilaments  at  the  edges  of  the  microtubule. 

7.4  Theoretical  Tubulin  Model  Details 

Figure  7.2a  is  a  schematic  of  our  tubulin  model.  The  model  consists  of  an  impen¬ 
etrable  rectangle  of  dimensions,  4.6  by  8.0  (unitless),  encasing  four  point  charges 
of  magnitude  1.  The  point  charges  are  arranged  such  that  the  model  protein  has 
zero  charge,  and  the  dipole  moment  is  directed  along  the  long  axis  of  the  protein 
(parallel  with  the  microtubule  axis).  The  charges  are  separated  by  vertical  and  lat¬ 
eral  distances,  s  and  t,  so  that  if  s  ^  0  or  t  ^  0,  moments  of  order  higher  than  the 
dipole  are  generated.  This  model  is  a  generalized  version  of  the  point  dipole  model 
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of  tubulin  used  by  Tuszynski  et  al.  [92] ,  and  it  reduces  to  a  point  dipole  model  when 


0  and  t  — >  0. 


7.5  Energetic  Considerations 

The  model  proteins  are  arranged  in  a  3D,  cylindrical,  microtubule-like  lattice  as 
shown  in  Figure  5.7b.  We  use  a  Coulombic  energy  equation, 

,  interprotein 

Hij 

0,  intraprotein 

and  sum  over  all  pairs  of  charges  (excluding  the  energy  of  pairs  of  charges 
within  the  same  protein)  to  obtain  the  energy  of  our  model  microtubule, 

4 zc  4 zc 

Urn  =  EE  Utj,  (7,1) 

i=  1  j=i+ 1 

where  i  is  a  parameter  that  runs  through  the  charges  along  the  length  of  the  mi¬ 
crotubule,  Z  is  the  number  of  proteins  in  each  protofilament,  j  is  a  parameter  that 
runs  through  the  charges  around  the  circumference  of  the  microtubule,  and  C  is 
the  number  of  protofilaments.  In  our  investigation,  the  parameters  varied  were  the 
protofilament  staggering,  6,  the  vertical  and  lateral  charge  separations,  s  and  t,  the 
number  of  fraying  proteins  in  each  frayed  protofilament  at  the  edge  of  the  micro¬ 
tubule,  N,  and  the  curvature  of  the  frayed  protofilaments,  A.  C  was  kept  fixed 
at  13,  while  Z  was  varied  from  10  to  100.  We  normalize  Utot  to  the  energy  of  13 
infinitely  separated  protofilaments  of  length  Z,  such  that 
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Fig.  7.2:  a)  A  schematic  of  our  tubulin  model.  Four  point  charges  are  separated  by  ver¬ 
tical  and  lateral  distances,  s  and  t,  and  encased  in  an  impenetrable  rectangle  of 
dimensions,  4.6  by  8.0  (unitless),  b)  The  model  proteins  are  arranged  in  a  3D, 
cylindrical,  microtubule-like  lattice. 
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where  E  is  the  normalized  energy,  and  U f,j  is  the  energy  of  an  isolated  protofilament 
of  length  Z.  We  also  calculate  an  energy  density,  which  is  the  energy  of  each  model 
protein  within  a  structure, 


U{z,c)  =  E  E  (7.3) 

fc=internal  charges  z=external  charges  kl 

where  the  k  is  a  parameter  that  runs  through  all  the  charges  within  a  single  pro¬ 
tein,  /  is  a  parameter  that  runs  through  all  charges  outside  of  that  protein,  z  is  an 
integral  position  index  of  the  protein  in  a  protofilament  (1  <  z  <  Z),  and  c  is  an 
integral  position  index  of  the  protofilament  within  the  microtubule  (1  <  c  <  C).  A 
normalized  energy  density, 


£(z,c) 


Utot{z,c) 

Mfil(z) 


(7.4) 


is  then  calculated,  where  Utot(z,  c)  is  the  energy  density  of  a  protein  at  position  (z,  c ) 
within  the  microtubule  of  length  Z,  and  Ufu(z)  is  the  energy  density  of  a  protein  at 
position  z  within  an  isolated  protofilament  of  the  same  length,  Z. 

Figure  7.3a  is  a  plot  of  E  vs.  5  for  a  point  dipole  model  of  tubulin  (s  — ^  0 
and  t  — ^  0),  in  a  model  microtubule  where  Z  =  30.  Note  that  the  energy  minimum 
is  at  5  =  4.3  nm,  unlike  physical  microtubules.  This  indicates  that  the  point  dipole 
model  does  not  reproduce  a  realistic  staggering.  Figure  7.4a  is  a  plot  of  S  within  the 
model  microtubule.  Note  £  >  0  everywhere  within  the  microtubule.  This  is  evidence 
this  simple  point  dipole  model  is  cannot  produce  a  stable  microtubule  structure. 

Figure  7.3b  is  a  plot  of  E  vs.  6  for  a  four  point  charge  model  of  tubulin  with 
s  =  7.2,  t  =  4.4,  and  Z  =  30.  Note  that  the  energy  minimum  is  at  5  =  0.93 
nm,  as  is  the  case  in  13-protohlament  microtubules.  Figure  7.4b  is  a  plot  of  £ 
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within  the  same  model  microtubule.  S  <  0  within  the  core  of  the  microtubule, 
indicating  the  formation  of  a  stable  microtubule  core.  The  sharp  rise  in  S  at  the 
ends  of  the  microtubule  is  evidence  that  the  microtubule  edges  are  unstable  and 
have  an  energetic  preference  to  separate  into  isolated  protofilaments.  The  four 
charge  model  of  tubulin  with  s  =  7.2  and  t  =  4.4  reproduces  both  the  0.93  nm 
staggering  of  real  microtubules,  and  the  stable  tube  structure.  The  instability  of 
the  model  microtubule  ends  also  suggests  that  the  model  may  reproduce  the  fraying 
protofilaments  seen  in  earlier  experimental  studies  [82], 


7.6  Conclusions 

In  this  chapter,  we  showed  that  if  the  interactions  of  all  particles  are  uniformly  al¬ 
tered,  complex  self-assembled  structures  can  be  stabilised.  As  an  example,  we  pro¬ 
posed  a  simple  point  charge  model  of  the  protein  tubulin,  and  showed  that,  without 
altering  the  dipole  moment,  the  addition  of  higher  order  multipole  moments  to  the 
model  can  energetically  stabilise  the  biologically  relevant  self-assembly  pattern.  We 
also  presented  a  physical  macroscopic  magnetic  analog  of  our  model,  and  show  that 
it  has  many  of  the  same  structural  characteristics  as  microtubules. 

7. 7  Outlook 


Given  that  a  simple  four  point  charge  model  can  stabilise  such  complex  structures 
as  microtubules,  it  is  likely  that  similar  models  could  be  used  to  gain  insight  into 
the  self-assembly  of  other  anisotropically  interacting  particles  and  macromolecules. 
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Fig.  7.3:  a)  A  plot  of  E  vs.  6  for  a  point  dipole  model  of  tubulin  in  a  3D  microtubule 
lattice  with  Z  =  30.  Note  that  the  energy  minimum  is  at  5  =  4.3nm,  unlike 
physical  microtubules,  b)  A  plot  of  E  vs.  5  for  a  four  point  charge  model  of 
tubulin  with  s  =  7.2,  t  =  4.4,  and  Z  =  30.  Note  that  the  energy  minimum  is  at 
5  =  0.93nm,  as  is  the  case  in  real  13-protofilament  microtubules. 
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Fig.  7.4:  a)  A  plot  of  £  within  the  model  microtubule.  Note  £  >  0  everywhere  within 
the  microtubule,  indicating  instability,  b)  A  plot  of  £  within  the  same  model 
microtubule.  £  <  0  within  the  core  of  the  microtubule,  indicating  the  formation 
of  a  stable  microtubule  core.  The  sharp  rise  in  £  at  the  ends  of  the  microtubule  is 
evidence  that  the  microtubule  edges  are  unstable  and  have  an  energetic  preference 
to  separate  into  isolated  protofilaments. 


97 


Developing  such  models  could  be  difficult,  however,  without  deriving  the  models 
directly  from  the  crystallographic  data  in  the  Protein  Data  Bank,  similar  to  the 
method  described  in  Chapter  6.  In  order  to  generalize  this  modeling  to  function 
for  any  arbitrary  protein,  methods  must  be  developed  to  effectively  and  simply 
model  the  hydrophobic  interactions  betwee  proteins.  Together  with  the  method  of 
modeling  electrostatics  described  in  Chapter  6,  a  method  of  developing  a  minimal 
hydrophobic  model  of  protein  interactions  could  lead  to  an  efficient  method  of  mod¬ 
eling  and  simulating  a  wide  array  of  proteins,  and  could  give  important  insight  into 
the  self-assembly  and  possibly  aid  in  drug  discovery. 
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APPENDIX 


A.  SEGREGATION  IN  A  MONOLAYER  OF  MAGNETIC 


SPHERES 

A.l  Energetic  Considerations 

We  assume  our  magnets  are  uniformly  magnetized  cylinders  of  axial  length,  /,  con¬ 
stant  diameter,  d,  with  the  cylinder  axis  coincident  with  the  z  axis  and  the  center 
of  mass  at  the  origin.  We  may  then  model  the  magnetostatics  of  the  magnet  using 
two  uniformly  charged  discs  of  opposite  sign,  but  equal  total  charge  magnitude  qtot , 
of  diameter  d  with  the  centers  of  each  disc  set  a  distance  l  apart,  and  with  the 
coincident  planes  of  each  of  the  discs  being  parallel.  The  magnetostatic  moment 
calculations  of  the  magnets  follow. 

Charge: 


Q  =  /  ,,  p(x)d3x  =  qtot  -  qtot  =  0 
Jail  space 


(A.l) 


Dipole: 


Pr  = 


/  all  space 


p(x)a;d3x  =  0 


(A.2) 


and 


(A.3) 


Py=[  p(x)|/rf3X  =  0, 

Jail  space 

by  symmetry. 

Pz=  p(x)^d3x  =  qtotl  -  qtot(-l)  =  2 qtotl  (A.4) 

Jail  space 

Quadrupole: 

Qxx  =  /„  p(x)(2a;2  -  y2  -  z2)d3x  =  0,  (A.5) 

J all  space 

because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z.  Simi¬ 
larly,  Qyy  QzZ  O’ 

Qxy  =  Qyx=  I  p(x)(3xy)d3x  =  0,  (A. 6) 

J  cill  Sp&CG 

because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 

Qxz  =  Qzx=[u  p(x)(3xz)d3x  =  0,  (A. 7) 

Jail  space 

because  the  position  term  is  odd  in  x  and  the  charge  term  is  even  in  x.  Simi¬ 
larly,  Qyz  =  Qzy  =  0,  by  symmetry. 

Thus,  all  quadrupole  moments  are  identically  zero  because  of  symmetry. 
Octopole: 

We  also  calculate  the  octopole  terms  in  the  multipole  moment  expansion  (equa¬ 
tion  2.1).  In  particular,  we  calculate  the  l  =  3  terms  in 
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Qlm  —  /  11  Yl*m(0,<l>)p(xyd3X, 

Jail  space 


(A.8) 


where  the  Y^n(9,(j>)  are  the  complex  conjugates  of  the  spherical  harmonics. 


qso  =  \/T^r  /  n  (22*3  ~  3-r2)p(x)rf3x 

V  lo7r  Jail  space 

=  f d/ 2  -  6r3)rfr 

7 rd2  V  167T  ir=0 

_  4/gtoi  /  7  l2d 2  _  6c?4 

c?2  V  167t  8  64 

Iqtot  I  7  3c?  \ 


(A.9) 


(A.  10) 


(A.  11) 


(A.  12) 
(A.  13) 


531  =  “IvS/all  space*1  “  "  rXxKix  =  °- 


(A.14) 


because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 


932  =  IV^T  /all  s  ace^2  ~~  ^  ~  2ia;^Mx)d3x 


space 


1  105  r 

4  V  27t  Jail  space 


-2  ixyz)p(x.)d3?£  =  0, 


(A.  15) 


(A. 16) 


because  the  position  term  is  odd  in  x  and  the  charge  term  is  even  in  x. 


1  /  35  r 

g33  =  --y  —  yaU  gpace(x3  +  ip3  -  3au/2  -  3ia;2?/)p(x)c?3x  =  0,  (A.17) 


because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 
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The  magnetostatic  moment  calculations  of  our  model,  as  described  in  2.3 
follow.  For  the  purposes  of  these  calculations,  the  coincident  plane  of  the  four 
charges  is  the  yz  plane,  with  the  center  of  mass  at  the  origin,  the  dipole  moment 
parallel  with  the  z  axis,  and  the  magnitude  of  the  charges,  q  =  to  keep  the  total 
dipole  moment  constant  and  equal  to  that  of  the  magnet. 

Charge: 


Di-pole: 


Q=  p(x)d3x  =  q  +  q  -  q  -  q  =  0  (A. 18) 

J all  space 


Px=  p(x)a;d3x  =  0 

Jail  space 


(A. 19) 


and 


Py  =  /  ii  p(x)r/d3x  =  0, 
Jail  space 


(A. 20) 


by  symmetry. 


P ,  = 


/  p(x)zd3x  =  qh  -  q(~h)  =  2 qh  =  2 qtotl  (A.21) 

Jail  space 


Quadrupole: 


Qxx  =  /  „  p(x)(2a;2  -  y2  -  z2)d3x  =  0, 

Jail  space 


(A. 22) 


because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z.  Simi¬ 


larly,  Qyy  Q  ZZ  O’ 
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Qxy  =  Qyx  =  I  „  p(x)(3xy)d3x  =  0,  (A. 23) 

J  8.11  Sp8CG 

because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 

Qxz  =  Qzx=lu  p(x.)(3xz)d3xi  =  0,  (A. 24) 

./all  space 

because  the  position  term  is  odd  in  x  and  the  charge  term  is  even  in  x.  Simi¬ 
larly,  Qyz  =  Qzy  =  0,  by  symmetry. 

Thus,  all  quadrupolc  moments  are  identically  zero  because  of  symmetry. 
Octopole: 

We  also  calculate  the  octopole  terms  in  the  multipole  moment  expansion  (equa¬ 
tion  2.1).  In  particular,  we  calculate  the  l  =  3  terms  in 


qim=f  Yl*rn(6,(j))p(x)rld3x,  (A.25) 

J  ail  space 

where  the  (j))  are  the  complex  conjugates  of  the  spherical  harmonics. 


.1  [  „  (2z3  -  3zx2  -  3 zy2)p{x)d3* 

1o7t  Jail  space 

1 - 

(A. 26) 

=  J  qi(2z3  3zx2  3zv2) 

'  71  all  charges 

(A. 27) 

=  9tc*hJ  7  (2 h2  3 w2) 

4  V  167tv 

(A. 28) 

(A. 29) 

1  /  21  r 

"31  =  “iVs/all  space** " iy)(5z  ~  "  )p(x)<*  X  =  (A30) 
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because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 


q32  =  7  y  X  (b(zx2  ~  ZV 2  “  Hxyz) 

’  77  all  charges 


1  / 105 
4V  2tT 


X!  <n(-zy2) 

all  charges 

qhw 2  / 105 

8  V  2tT’ 


(A.31) 
(A. 32) 
(A. 33) 


1  /  35  r 

^33  =  Ty  /  „  (a:3  +  *l/3  -  3au/2  -  3 ix2y)p(x)d3x  =  0,  (A.34) 

4  V  47 r  Jail  space 

because  the  position  term  is  even  in  z  and  the  charge  term  is  odd  in  z. 

In  matching  the  model  to  the  experiment,  all  moments  below  octopolc  are 
equivalent  by  design.  For  the  octopolc  terms, 


(A. 35) 
(A. 36) 

Thus,  h  —  0  or  w  —  0  for  the  model. 

Also, 


while  <732  = 


<?32  =  0  for  the  experiment, 
qhw 2 


'105 

2 V 


for  the  model. 


<?30  — 

while  qso 


lq, 


tot 


167T 


0 h 2 


Qtoth  /  7 


167T 


for  the  experiment, 
3 w2)  for  the  model. 


Thus, 


(A. 37) 
(A. 38) 


105 


<?30  — 


So,  if 


and  for 


iQtot 

~Y 


3  “ 


l  >  w  =  0 

<  h  =  0 

l  = h  =  w  =  0 

h  =  °,  w  =  v/f  -  ¥ 

w  =  0,  h  =  yj l2  - 
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a. 


B.  PATTERN  FORMATION  IN  A  MONOLAYER  OF 


MAGNETIC  SPHERES 

B.l  Energetic  Considerations 

Figure  B.l  is  a  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely 
separated  chains  of  length  y/N)  of  several  different  patterns  of  short  (h  =  0.001, 
w  =  0.25)  model  particles  versus  y/N .  All  patterns  except  the  ring  pattern  are 
rectangular  lattices  of  dimensions  N  particles  by  N  particles.  Note  that  the  mi¬ 
crovortex  and  macrovortex  energies  are  nearly  identical  and  are  the  energetically 
preferred  patterns  for  \/N  >  3.  Also  note  that  the  HCP  ferromagnetic  pattern 
decreases  significantly  relative  to  all  other  patterns,  over  the  calculated  range. 

Figure  B.2  is  a  plot  of  the  energies  (normalized  to  the  energy  of  y/N  infinitely 
separated  chains  of  length  y/N)  of  several  different  patterns  of  weak  (h  =  0.001,w  = 
0.31)  model  particles  versus  y/N.  All  patterns  except  the  ring  pattern  are  rectangu¬ 
lar  lattices  of  dimensions  N  particles  by  N  particles.  Note  that  the  microvortex  is 
the  energetically  preferred  pattern  for  y/N  >  3.  Also  note  that  the  HCP  ferromag¬ 
netic  pattern  decreases  significantly  relative  to  all  other  patterns,  over  the  calculated 


range. 
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Fig.  B.l:  A  plot  of  the  energies  (normalized  to  the  energy  of  \fN  infinitely  separated  chains 
of  length  \/N)  of  several  different  patterns  of  h  =  0.001,  w  =  0.25  particles  versus 
sfN .  All  patterns  except  the  ring  pattern  are  rectangular  lattices  of  dimensions 
N  particles  by  N  particles.  Note  that  the  nricrovortex  and  macrovortex  energies 
are  nearly  identical  and  are  the  energetically  preferred  patterns  for  y/N  >  3. 
Also  note  that  the  HCP  ferromagnetic  pattern  decreases  significantly  relative  to 
all  other  patterns,  over  the  calculated  range. 
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Fig.  B.2:  A  plot  of  the  energies  (normalized  to  the  energy  of  \/fV  infinitely  separated 
chains  of  length  \/N)  of  several  different  patterns  of  h  =  0.001, re  =  0.31  parti¬ 
cles  versus  y/~N.  All  patterns  except  the  ring  pattern  are  rectangular  lattices  of 
dimensions  N  particles  by  N  particles.  Note  that  the  nricrovortex  is  the  ener¬ 
getically  preferred  pattern  for  y/~N  >  3.  Also  note  that  the  HCP  ferromagnetic 
pattern  decreases  significantly  relative  to  all  other  patterns,  over  the  calculated 
range. 
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C.  PATTERN  FORMATION  IN  A  MONOLAYER  OF 
MAGNETIC  SPHERES 

C.l  Energetic  Considerations 

Below  are  the  equations  used  in  the  magnetostatic  energy  calculating  algorithm. 
The  calculations  are  performed  in  Mathematica.  In  general,  the  calculations  follow 
the  following  steps: 

1.  Determine  charge  positions,  r  in  the  pattern  being  investigated. 

2.  Determine  distances  between  any  two  charges,  R ,  in  the  pattern. 

3.  Determine  charge  values,  Q. 

4.  Calculate  the  magnetostatic  energy,  U  of  the  entire  pattern.  Note:  this  step 
requires  careful  attention  to  summation  limits  so  as  not  to  double  count  or 
include  particle  self-energies. 

The  calculation  of  a  linear  chain  begins  with  the  determination  of 
all  charge  positions, 

I  (0,  ^  —  §),  rn  even 
rune(d,  s,m)  =  < 

1  (Q;  irn-l)d  _)_  m  odd 


where  m  is  a  dummy  variable  unique  to  each  charge.  The  distances  between  all 
pairs  of  charges  is  then  simply 


Rune{d i  s,  m,  n)  =  rUne(d ,  s,  m)  -  rnne(d,  s,  n) 


(C.l) 


where  n  is  another  dummy  variable  unique  to  each  charge.  A  charge  value,  Qune{q ,  d,  s,  m) 
is  then  assigned  to  each  charge  such  that 


Qune(q,d,s,m )  = 

The  energy  of  the  linear  chain  is  then 


-  — ,  m  even 

.s  * 


m  odd 


2x _ 3  2X _ 1  I  (1  fi?n,n  —  l) f-iQ  Qiine(1hd',s,7n^Qiine(q,d,s,n)  even 

TT  \  \  '  \  '  I  4tt  \RUne,(d,s,m,n)\  '  ‘ 

Uune(X,  s,  d,  q)  =  2^  2^  \ 

m= 0  n=m+ 1  |  n0  Qune(q,d,s,m)Qiine(q,d,s,n)  m 


4tt  \Rune(d,s,m,n)\  ’ 

where  5TO)„-i  is  the  Kronecker  delta,  and  X  is  the  number  of  particles  in  the  chain. 

Similarly,  the  positions  of  all  charges  in  an  HCP  ferromagnetic  pat¬ 
tern  can  be  calculated, 


rhex(d,s,m,n)  = 


2  ’ 

2  2 

\/3  dn 

(m—l)d 

2  ’ 

2 

y/3dn 

(m+l)d 

2  ’ 

2 

yj  3dn 

0  1 

md  .  s 
0  i  0 

m  even,  n  even 


m  odd,  n  odd 

where  (m,  n)  is  a  set  of  dummy  variables  unique  to  each  charge.  The  distances 
between  all  pairs  of  charges  is  then  simply 


Rhex(d,  s,  m,  n,  /,  k)  =  rhex(d,  s,  m,  n)  -  rhex(d,  s,  /,  k) 


(C.2) 
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where  (/,  k)  is  another  set  of  dummy  variables  unique  to  each  charge.  A  charge 
value,  Qhex(q,d,  s,m)  is  then  assigned  to  each  charge  such  that 


Qhex(q,d,s,m ) 


qd 
s  ’ 


m  even 


The  energy  of  the  rectangular  HCP  ferromagnetic  lattice  of  dimensions  X  particles 
x  Y  particles  is  then 


Uhex i^X .Y^s^d^q)  YUune(X j  5,  d,  q)  T 

IM)  Qhex(q,  d,  s,  m)Qhex(q,  d,  s,  l) 

io  h  kY+\  4n  I  Rhex{d,  s,  m,  n,  l,  k )  | 


(C.3) 


Also,  the  positions  of  all  charges  in  a  square-packed  antiferromag¬ 
netic  pattern  can  be  calculated, 

rantsq(d,  s,  m,  n)  = 


( dn , 

md  s  \ 

2  2  /  ’ 

m  even 

( dn , 

(m-l)d  |  s\ 

2  "T"  2 /’ 

m  odd 

where  (m,  n)  is  a  set  of  dummy  variables  unique  to  each  charge.  The  distances 
between  all  pairs  of  charges  is  then  simply 


kdantsq(di  A  A  Y’antsq(d')  S,  /7l,  n)  ^antsq(^d,  S,  Z,  /c) 


(C.4) 


where  (/,  k )  is  another  set  of  dummy  variables  unique  to  each  charge.  A  charge 
value,  Qantsq(q,d,s,m,n)  is  then  assigned  to  each  charge  such  that 

— ,  m  even,  n  even 

—  — ,  m  odd,  n  even 

—  — ,  m  even,  n  odd 
— ,  m  odd,  n  odd 

.s  * 


Qantsq(q,d,s,m,n)  K 
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The  energy  of  the  rectangular  square-packed  antiferromagnetic  lattice  of  dimensions 
X  particles  x  Y  particles  is  then 


Uantsq(X,  V,  s,  d,  (/)  1"  Unne(X,  s,  d,  (/)  T 

^ ^  ^ ^  ^ ^  /7q  Qant.sq^Qi  d,  S,  U / .  Tl)Q antsqiSl i  d-  S,  /,  /c) 


m=0  n=0  (=0  fc=n+l 


47T 


Rantsq(d,  S ,  / / / .  '77,  t,  /c) 


(C.5) 


The  positions  of  all  charges  in  an  HCP  macrovortex  pattern  can  be 
calculated,  using  some  additional  equations  to  account  for  the  orienta¬ 


tions  of  the  particles, 


I'Xmacro(d)  Tl ) 


v/3d 


n 


- — m  even,  n  even 


^ Y macro  (d ■,  TH)  Tl^)  \ 


md 
2  ’ 


Qmacro(d,m,n )  =  tan 


m  odd, n  even 

—  uy,  m  even,  n  odd 

(m~1-)d,  m  odd,  n  odd 
_i(X  macro 


Y  macro ' 


Finally,  the  positions  of  the  charges  are  determined  by 


(C.6) 


(C.7) 


^macro^d,  S,  777.,  77  ) 

( Y^dn  _|_  ^"1+1)d+|sm(emac7.0(d,?7i,?i))),  ?77  even,  n  even 

(  x/itdra  _j_  |cos^g 

macro  (d,m,n)),  +  ^sm(7r+0  macro  (d,m,n))),  777  odd,  n  even 

( — 7T^  +  ^cos{n+9rnacro{d,m,n)),n^L  +  ^sin{6macro(d,m,n))),  777  even,  Tl  odd 

(  \ZHdn  _|_  | COs(ernacro(d,m,n)),  2^d  +  %sin(Tr+9macro(d,m,n))),  777  odd,  '71  odd 


(C.8) 
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where  (m,  n)  is  a  set  of  dummy  variables  unique  to  each  charge.  The  distances 
between  all  pairs  of  charges  is  then  simply 

Rmacroidi  ^b  ^ macro(di  ^b  ^0  ^"macroidi  S ,  Z,  /c)  (C.9) 

where  (/,  k)  is  another  set  of  dummy  variables  unique  to  each  charge.  A  charge 
value,  Qmacroiq,  d,  s,  m)  is  then  assigned  to  each  charge  such  that 

{—  y,  m  even 
m  odd 
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The  energy  of  the  HCP  macrovortex  lattice  of  dimensions  X  particles  x  Y 


particles  is  then 


X  X 

Y- 1  T  T 

E  E  E 


Y—l 

^  MO  Qmacro(Qid,S,7n^Qmacro{Q,d,S,l) 


- ,  X,Y  even ,  n^k 


T=-y  n=_x=l  k=n 


I y  47T  |  Rmacro  (d,S,m,n,l,k^  | 

y—i 

EMO  Q macro  {Qid,S,m,')Qrnacro(Qid,S,l) 

- ,  X,Y  even ,  m  odd ,  n=/c 

4-7T  |i^macro(d,S,m,n,/,fe)| 

y—i 

^  •>  MO  Qmacro(Qid,S,m,')Qrnacro(q,d,S,l) 


) 

-j  X,Y,m  even ,  n=/c 


y-i 


X  X 

2  2 


E  E  E 


i=-y+l  k=n 


(  l=m+ 2  4?r  l^macro(rf,S,Wl,^,fe)| 

f  y—i 

J  \  ■>  MO  Qmacro{Qjd,S,Tn)Qmacro{Q^d,S,l^ 

I  y - ,  X  even,  Y  odd,  n^k 

l _ _ y_|_l  4-7T  |  Rmacro  (d,S,?7i,7l,Z,/c)  \ 

Y 

^  MO  Qmacro(Qid,S,Tn^Qrnacro(Qid,S,l') 

» - ,  X,y  even,  Y,m  odd,  n- 


E 

Z=m+1 


4-7T  | -Rmacro  (d, S, 777, 77, Z,fc)| 


y 

EMO  Q macro  ( Q:d,S,m')Qrnacro  (,Q)d,S,l') 

- ?  X,Y,m  even,  Y  odd,  n- 

Astt  IP  Ml 


X-l  X-l 

y—l  2  2 

E  E  E 

m=—Yn—  X~1  k=n 


^  2  4?r  | Rmacro  (d,s,m,n,l,h)  | 

f  y—i 

EMO  Qmacro{Qid,S,7n')Qmacro(Q,d,S,l ) 

- ,  y  even,  X  odd,  n^k 

£ _ _ y  47T  | Rmacro  (d,S,Tn,n,l,k')  | 

y—i 

^  •>  MO  Qmacro(Qid,S,m,')Qrnacro(Qid,S,l') 
l=m+l  477  \Rmacro(d,s,m,n,l,k)\ 

Y—l 

^  MO  Qmacro(Qid,S,m,')Qmacro(Qid,S,l') 

^  4?r  |  Rmacro  (d, S, 777, 77, Z,fc)| 


,  y  even,  X,m  odd,  n=k 

,  y,m  even,  X  odd,  n=k 


Y 


x-i  x-i 

2  2 


E  E  E 

m=— y+1  n— _  2£n!  k=n 


'  y 

EMO  Q  macro  ( Q^d,S,7n^Qrnacro{q^d,S,l ) 

- ,  X,y  odd,  n//c 

^ y_|_^  47T  [Rmacro  (d,S, 777, 77, Z,/c)  | 

y 

EMO  Q  macro  ( Q:d,s,m')Qrnacro(^q,d,s,l ) 

- ,  X,Y,m  odd,  n=k 

l=m+ 1  4-7r  l-Rmocro(d,s,m,n,Z,fe)| 

y 

^  MO  Qmacro(si,)d,S,ni)Qrriacro(si’>d,S,V) 

^  ^_m_|_2  4?r  |Rmacro(d,S,777,77,Z,/c)| 


d) 

— ,  m  even,  X,Y  odd,  n- 


(C.10) 
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The  positions  of  all  charges  in  a  square-packed  microvortex  pattern 


can  be  calculated, 

(< in  —  yp,  y  ~  m,  n,  y  even 

(dn+^^-f  -  m,  n  even,  f  odd 

( dn  +  +  pp),  n,  even,  m  odd 

(dn  —  pp,  (m~1)d  +  pp),  n  even,  m,  odd 

rmicro(d,s,m,n)  = 

(dn— pp,f^  +  pp),  m,  f  even,  n  odd 

(dn  +  pp,  'y  +  pp),  m  even,  n,  y  odd 

(dn  +  pp,  -  pp),  pp  even,  m,  n  odd 

(dn  -  pp,  -  pp),  m,  n,  pp  odd 

where  (m,  n)  is  a  set  of  dummy  variables  unique  to  each  charge.  The  distances 
between  all  pairs  of  charges  is  then  simply 

ddmicro(d')  5,  UT,  U,  /,  /c)  V micro(d-!  5,  m,  U.)  r micro(di  5,  /,  /t)  (C.ll) 

where  (/,  k)  is  another  set  of  dummy  variables  unique  to  each  charge.  A  charge 
value,  Qmicro(q,d,s,m )  is  then  assigned  to  each  charge  such  that 

{— y,  m  even 
y,  m  odd 
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The  energy  of  the  square-packed  microvortex  lattice  of  dimensions  X  particles 


x  Y  particles  is  then 


Umicro  ( X ,  T,  S ,  rf,  (/) 

'  Y- 1 

EMO  Q  micro  (q,d,S,m')Q micro  (.Q)d,S,V) 

- ,  X,y  even,  n^k 

j y  4tt  | Rmicro {d,s,m,n,l,k') | 

y_l  f-  y_i 

Ey  y  y  MO  Qmicro(Qid,s,m)Qrnicro(q,d,s,l) 

y  /  i  / - ">  X,Y  even ,  m  odd,  n=k 

m=-Y  n=_X-2  k=n  l=m+ 1  4w  \Rmicro(d,s,m,n,l,k)\ 

2  V— 1 

EMO  Q micro  (Sl^d,S,TTl)Q'fYiic<ro  (yQ)d,S,V) 

- ,  X,Y,m  even,  n=k 

l=m+ 2  4?r  |-Rmicro(^,5,m,n,/,fc)| 

y 

EMO  Q micro  (Qid,s,Tn')Qrnj/Cro  i^q,d,s ,/) 

- ,  X  even,  Y  odd,  n^k 

l=_yjrl  4-7T  |i?micro(rf,s,m,n,/,fc)| 

v  y  -v  MO  QmicroiQid^jm^Qmicro^qjdjS, l) 

Y  y  y - ,  X,Y  even,  Y,m  odd,  n=k 

vi=-Y+ln=_x=lk=n  l=m+ 1  4w  \Rmicro(d,s,m,n,l,k)\ 

2  Y 

EMO  Qmicro(Qid,S,Tn')Qrnicro(q,d,S,l ) 

- ,  X,Y,m  even,  Y  odd,  n=k 

4^"  | Rmicro (d,s,yn,n,l,k') | 

'  y_l 

EMO  Q micro  (,Qid,s,m')Qrnicro  (q,d,s,V) 

- ?  y  even,  X  odd,  n^k 

i y  477  |  Rmicro  {d,S,Tn,n,l,k')  \ 

X  —  1  X  —  1  L  J 

Ey  "v  y  -v  y  >  MO  Qmicro{Q,d,S,m)Qrnicro(q,d,S,l) 

y  y  \  y - ,  y  even,  X,m  odd,  n=k 

rn=—Yr)—_K=l  k=n  I  /=m+l  47r  l^tero(d,«,m)n)i,fc)j 


EMO  Q micro  (,Q id, S,Tf\)Q micro  (^Q^d,S ,Z) 

- ,  y,m  even,  X  odd,  n=k 

l=m+ 2  4?r  l^m<cro(d,«,m,n,l,fe)| 

y 

EMO  Q micro  (q,d,S,m')Q micro  (.Q)d,S,V) 

- ,  X,y  odd,  n/A: 

i _ _ v_i_i  47t  \Rmicro(d,s,wn,n,l,k>)\ 

X-l  X-l  L~  1  _r± 

Y  ~2~  ~2~  Y 

Ev — "v  v—v  v — >  MO  Qmicro(Q,d,s,m)Qrnicro(q,d,s,l ) 

y  y  y - ,  X,Y,m  odd,  n=k 

i=-y+l  n-  x-l  k=n  l=m+ 1  4w  I Rmicro(.d,s,m,n,l,kj\ 

2  Y 

EMO  Q micro  ( Qid,s,m)Qrnicro  ( q,d,s,l ) 

- ,  m  even,  X,Y  odd,  n=k 

l=m+ 2  4?r  l#micro(d,s,m,ri,Z,/c)| 


(C.12) 


The  positions  of  all  charges  in  a  ring  pattern  can  be  calculated,  using 
an  additional  equation  to  account  for  the  orientations  of  the  particles, 


^ring^X^ 


n  odd 
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where  X  is  the  number  of  particles  in  the  ring.  Finally,  the  positions  of  the  charges 
are  determined  by 


T- ring(Xj  71,  d,  s) 


( _ i _ 

'2  s»(f) 


V2  sm(§) 


Sin(6ring  (-V 2  ring  )  > 

sin(9ring  (X,n))—  §  cos(6rirlg  ( X,n ) ) , 


2sin(  )  COS(^ring  (-V,7l))  2  ring  (  ) , 

2cin(  ^^i^ring  (-V,7l) )  “t"  2  ring  } 


n  even 

77  odd 


(C.13) 


where  n  is  a  dummy  variable  unique  to  each  charge.  The  distances  between  all  pairs 
of  charges  is  then  simply 


Faring  (d,  S,  777  ,  77,  -A)  7"rmg(-A,  777,  d,  s)  T"ring  (A",  ^7,  d,  s)  (C.14) 


where  777  is  another  dummy  variable  unique  to  each  charge.  A  charge  value,  Qring(q,  d,  s,  77) 
is  then  assigned  to  each  charge  such  that 


Qring^Qi  d,  S,  77) 


qd 

T  ’ 


77  even 


—  77  odd 

S  ’ 
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The  energy  of  the  ring  is  then 


U ringi^X ,  S,  rf,  q'j 


Xq2d2ido 
Ait  s3 


2X-2  2X-1 


+  E  E 

m= 0  n=m-\- 1 


A^O  Qring^Qi  d,  5,  'Wl'jQringi^Qi  d,  5,  Ti) 

47 r  I  Rring(d,  S,m,Tl,X)  \ 


(C.15) 
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